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1.  SUMMARY 


1.1  OVERVIEW 

This  project  addresses  the  problem  of  incorporating  seismic  wave  attenuation  into 
numerical  modeling  codes  which  use  time  stepping  to  integrate  the  equations  of  motion. 
Such  codes  include  time  domain  finite  difference  and  finite  element  methods.  CTBT 
verification  requires  an  understanding  of  the  propagation  of  regional  seismic  phases  over 
complex  regional  paths  which  cut  across  major  structural  boundaries.  The  computation  of 
synthetic  seismograms  by  finite  difference  methods,  2D  and  3D,  plays  an  important  role  in 
developing  such  understanding.  In  order  for  ground  motion  estimates  to  be  realistic, 
however,  the  models  must  account  not  only  for  regional  elastic  structure  of  the  path,  but 
also  for  anelastic  losses;  in  the  case  of  high-amplitude,  near  source  parts  of  the  propagation 
path,  nonlinear  losses  influence  the  signal  as  well,  and  a  comprehensive  numerical  model 
should  be  able  to  simulate  realistic  nonlinear  attenuation. 

We  report  work  done  under  Grant  No.  F49620-94- 1-0205  during  the  period  March 
1994  through  March  1997.  This  work  has  been  done  in  coordination  with  a  parallel  effort 
led  by  J.  Bernard  Minster  at  the  University  of  California,  San  Diego  (Grant  F49620-94-1- 
0204).  Section  2  addresses  the  computational  problem  of  modeling  realistic  linear, 
anelastic  attenuation  in  2D  and  3D.  Heretofore,  3D  finite  difference  simulations  have 
generally  neglected  attenuation,  due  in  large  part  to  the  onerous  storage  requirements 
entailed  by  realistic  seismic  Q  models.  We  describe  a  novel  solution  to  this  problem  which 
we  call  memory-variable  coarse-graining.  The  coarse-graining  method  leads  to  an  order  of 
magnitude  reduction  in  computer  storage  requirements  for  anelastic  memory  variables. 

Section  3  develops  and  demonstrates  a  method  for  modeling  nonlinear  wave 
propagation  in  rock  under  conditions  of  intermediate  strain,  defined  as  strain  levels  below 
the  threshold  for  failure  and  damage,  but  above  the  threshold  for  onset  of  nonlinear 
response.  We  develop  a  model  which  reproduces  the  stress-strain  behavior  of  rock  in  this 
strain  regime,  as  measured  in  quasi-static  laboratory  tests,  and  show  that  the  resulting 
model  is  consistent  with  key  features  of  laboratory  wave  propagation  measurements. 

1.2  SUMMARY  OF  SECTION  2  -  ANELASTIC  COARSE-GRAINING 

Improvements  in  computing  speed  have  progressively  increased  the  usable 
bandwidth  of  seismic  wavefield  simulations  computed  with  time-stepped  numerical 
schemes  (e.g.,  finite  difference,  finite  element,  pseudospectral).  As  computational 
bandwidth  increases,  anelastic  losses  become  increasingly  significant  for  some  important 
applications  such  as  earthquake  ground  motion  modeling,  whole  earth  seismogram 


1 


simulation,  and  exploration  seismic  profile  modeling,  and  these  losses  need  to  be  included 
in  the  simulations.  As  bandwidth  increases,  however,  the  memory  variables  necessary  to 
incorporate  realistic  anelastic  losses  account  for  an  increasing  proportion  of  total 
computational  storage  requirements,  a  consequence  of  the  broad  relaxation  spectrum  of 
typical  earth  materials.  To  reduce  these  storage  requirements,  we  introduce  a  new  method 
in  which  the  memory  variables  are  coarse-grained,  i.e.,  redistributed  in  such  a  way  that 
only  a  single  relaxation  time  is  represented  at  each  node  point  (and  therefore  a  single 
memory  variable  per  stress  component  is  required).  Guided  by  a  perturbation  analysis,  we 
effect  this  redistribution  in  such  a  way  that  spatial  variability  of  this  single  relaxation  time 
simulates  the  full  relaxation  spectrum.  Such  coarse-graining  reduces  memory-variable 
storage  requirements  by  a  factor  of  8  for  3D  problems,  or  a  factor  of  4  for  2D  problems. 

In  fourth-order  finite  difference  computations  for  the  3D  acoustic  wave  equation, 
the  method  simulates  frequency-independent  Q  within  a  3%  tolerance  over  2  decades  in 
frequency,  and  is  highly  accurate  and  free  of  artifacts  over  the  entire  usable  bandwidth  of 
the  underlying  finite  difference  scheme.  These  results  should  also  hold  for  the 
elastodynamic  equations.  The  method  is  readily  generalized  to  approximate  speciric 
frequency-dependent  Q  models  such  as  power  laws,  or  to  further  reduce  memory 
requirements.  In  its  present  implementation,  the  main  limitation  of  the  method  is  that  it 
generates  artifacts  at  wavelengths  equal  to  4  grid  cell  dimensions  and  shorter,  which  may, 
in  some  limited  circumstances,  overlap  the  usable  bandwidth  of  very  high-order  finite 
difference  and/or  pseudospectral  schemes. 

1.3  SUMMARY  OF  PART  II  -  NONLINEAR  WAVE  PROPAGATION 

We  develop  a  method  for  modeling  nonlinear  wave  propagation  in  rock  at 
intermediate  strain  levels,  i.e.,  strain  levels  great  enough  that  nonlinearity  cannot  be 
neglected,  but  low  enough  that  the  rock  does  not  incur  macroscopic  damage.  The 
constitutive  model  is  formulated  using  a  singular-kernel  endochronic  formalism,  and  this 
formulation  is  shown  to  satisfy  a  number  of  general  observational  constraints,  including 
producing  a  power  law  dependence  of  attenuation  (Q-1)  on  strain  amplitude.  Once  the 
elastic  modulus  is  determined,  and  a  second  parameter  fixed  to  give  linear  in  the  strain 
amplitude,  the  model  has  2  remaining  free  parameters.  One  of  these  represents  cubic 
anharmonicity,  and  we  set  in  to  agree  with  laboratory  observations  of  harmonic  distortion. 
The  other  parameter  controls  the  amount  of  hysteresis,  and  it  is  set  to  approximate  stress- 
strain  curves  measured  in  laboratory  uniaxial  stress  experiments  on  Berea  sandstone  by 
Boitnott  and  Haupt.  The  constitutive  equations,  though  fundamentally  nonlinear  and  rate- 
independent,  have  a  superficial,  formal  resemblance  to  viscoelasticity,  which  we  exploit  to 
produce  an  efficient,  stable  numerical  algorithm.  We  solve  ID  wave  propagation  problems 
for  this  constitutive  model  using  both  finite  difference  and  pseudospectral  methods.  These 
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methods  are  shown  to  reproduce,  to  high  precision,  analytical  results  for  quasi-harmonic 
wave  propagation  in  a  nonlinearly  elastic  medium. 

Application  of  the  Berea  sandstone  model  to  quasi-harmonic  wave  propagation 
shows  several  departures  from  results  obtained  with  nonlinear  elasticity.  The  Berea  model 
shows  more  rapid  decay  with  distance  of  the  fundamental  frequency  component,  due  to 
nonlinear,  amplitude-dependent  attenuation,  than  does  nonlinear  elasticity.  The  Berea 
model  also  shows  enhanced  excitation  of  the  order  3  harmonic,  in  agreement  with 
laboratory  observations.  In  addition,  the  growth  with  propagation  distance  of  the 
harmonics  of  the  source  excitation  shows  a  saturation,  relative  to  the  nonlinear  elasticity 
results.  This  behavior  reflects  the  competing  effects  of  amplitude  growth  via  energy 
transfer  from  the  source  frequency,  and  energy  dissipation  due  to  hysteresis,  the 
dissipation  increasing  as  the  harmonic  amplitude  grows.  In  additional  numerical 
experiments,  we  find  that  a  two-frequency  source  function  generates  harmonics  with 
frequencies  which  can  be  expressed  as  linear  combinations  of  integer  multiples  of  the  two 
source  frequencies,  in  agreement  with  published  laboratory  results  for  other  solids. 
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2.  EFFICIENT  SIMULATION  OF  CONSTANT  Q  USING 
COARSE-GRAINED  MEMORY  VARIABLES 

Steven  M.  Day 


2.1  ABSTRACT 

Improvements  in  computing  speed  have  progressively  increased  the  usable 
bandwidth  of  seismic  wavefield  simulations  computed  with  time-stepped  numerical 
schemes  (e.g.,  finite  difference,  finite  element,  pseudospectral).  As  computational 
bandwidth  increases,  anelastic  losses  become  increasingly  significant  for  some  important 
applications  such  as  earthquake  ground  motion  modeling,  whole  earth  seismogram 
simulation,  and  exploration  seismic  profile  modeling,  and  these  losses  need  to  be  included 
in  the  simulations.  As  bandwidth  increases,  however,  the  memory  variables  necessary  to 
incorporate  realistic  anelastic  losses  account  for  an  increasing  proportion  of  total 
computational  storage  requirements,  a  consequence  of  the  broad  relaxation  spectrum  of 
typical  earth  materials.  To  reduce  these  storage  requirements,  we  introduce  a  new  method 
in  which  the  memory  variables  are  coarse-grained,  i.e.,  redistributed  in  such  a  way  that 
only  a  single  relaxation  time  is  represented  at  each  node  point  (and  therefore  a  single 
memory  variable  per  stress  component  is  required).  Guided  by  a  perturbation  analysis,  we 
effect  this  redistribution  in  such  a  way  that  spatial  variability  of  this  single  relaxation  time 
simulates  the  full  relaxation  spectrum.  Such  coarse-graining  reduces  memory-variable 
storage  requirements  by  a  factor  of  8  for  3D  problems,  or  a  factor  of  4  for  2D  problems. 

In  fourth-order  finite  difference  computations  for  the  3D  acoustic  wave  equation, 
the  method  simulates  frequency-independent  Q  within  a  3%  tolerance  over  2  decades  in 
frequency,  and  is  highly  accurate  and  free  of  artifacts  over  the  entire  usable  bandwidth  of 
the  underlying  finite  difference  scheme.  These  results  should  also  hold  for  the 
elastodynamic  equations.  The  method  is  readily  generalized  to  approximate  specific 
frequency-dependent  Q  models  such  as  power  laws,  or  to  further  reduce  memory 
requirements.  In  its  present  implementation,  the  main  limitation  of  the  method  is  that  it 
generates  artifacts  at  wavelengths  equal  to  4  grid  cell  dimensions  and  shorter,  which  may, 
in  some  limited  circumstances,  overlap  the  usable  bandwidth  of  very  high-order  finite 
difference  and/or  pseudospectral  schemes. 
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2.2  INTRODUCTION 


As  computing  has  become  progressively  cheaper,  it  has  become  possible  to  use 
discrete  numerical  methods  such  as  the  finite  difference  and  finite  element  methods  to 
synthesize  seismic  wavefields  over  spatial  domains  which  are  relatively  large  compared 
with  the  minimum  numerically  resolvable  wavelength.  When  a  seismic  pulse  propagates  a 
distance  many  times  its  dominant  wavelength,  anelastic  attenuation  cannot  generally  be 
ignored.  Even  three  dimensional  elastodynamic  problems  can  now  be  solved  over  domains 
sufficiently  large  that  attenuation  should  not  be  neglected,  and  the  problem  of  modeling 
attenuation  efficiently  and  accurately  will  become  increasingly  important  as  computational 
capabilities  advance.  Accurate  treatment  of  anelastic  attenuation  has  already  become  a 
priority  in  some  important  practical  applications  in  3D,  including,  for  example,  the 
modeling  of  earthquake  strong  ground  motion  in  sedimentary  basins  (e.g.,  Frankel  and 
Vidale,  1992;  Graves,  1993;  Yomogida  and  Etgen,  1993;  Olsen  et  al.,  1995a),  synthesis  of 
the  elastic  response  of  whole  earth  models  (e.g.,  Yoon  and  McMechan,  1995;  Igel  and 
Weber,  1995)  and  the  simulation  of  seismic  reflection  profiles  (e.g.,  Carcione  et  al.  1992). 
Applications  such  as  these  now  frequently  entail  three-dimensional  computations  over 
domains  with  dimensions  roughly  two  orders  of  magnitude  larger  than  the  minimum 
wavelength  (e.g.,  Olsen  et  al.,  1995b). 

The  most  general  anelastic  law  represents  the  stress  at  (x,0  as  a  convolution  integral 
over  the  complete  strain  history  at  x,  a  form  of  the  stress-strain  relation  which  is  not 
tractable  for  numerical  calculations,  since  it  imposes  enormous  storage  and  computational 
requirements.  Day  and  Minster  (1984)  developed  a  framework  for  incorporating  anelastic 
attenuation  laws  into  time-stepped  numerical  methods  such  as  the  finite  difference  method 
by  approximating  them  using  internal,  or  "memory",  variables.  This  approach  can  be 
viewed  in  the  time  domain  as  the  replacement  of  the  convolution  operator  by  a  low-order 
differential  operator.  Alternatively,  it  can  be  viewed  in  the  Laplace  transform  domain  as  the 
approximation  of  an  algebraic  transfer  function  by  a  pole- zero  representation.  Each  pole 
then  corresponds  to  a  memory  variable  which  evolves  according  to  a  first-order  differential 
equation  analogous  to  that  governing  a  standard  linear  solid.  The  latter  can  be  time-stepped 
to  compute  the  evolution  of  the  memory  variable(s)  along  with  the  other  field  variables  such 
as  velocity  and  stress.  This  framework  has  been  refined  and  improved  by  a  number  of 
investigators,  who  have  used  a  variety  of  methods  to  determine  efficient  pole-zero 
representations  of  the  stress-strain  relationship  (e.g.,  Emmerich  and  Korn,  1987;  Carcione 
et  al.,  1988,  Witte  and  Richards,  1990;  Krebs  and  Quiroga-Goode,  1994;  Blanch  et  al., 
1995). 

While  the  memory  variable  framework  is  effective,  it  does  impose  stringent 
demands  for  computer  memory.  To  approximate  an  anelastic  quality  factor  (Q)  which  is 
approximately  constant  over  a  specified  frequency  band  requires  2  to  3  memory  variables 
per  decade  of  bandwidth  per  stress  component  per  computational  unit  cell,  even  under  the 
most  highly  optimized  memory  variable  formulations  (e.g.,  Robertsson  et  al.,  1994; 
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Blanch  et  al.,  1995).  This  requirement  becomes  particularly  onerous  in  3D,  where  there 
are  6  independent  stress  components  per  computational  unit  cell.  The  typical  3D  elastic 
(non-attenuative)  staggered  grid  finite  difference  method  requires  9  words  of  storage  per 
unit  cell  for  the  field  variables,  i.e.,  the  velocity  and  stress  components.  Anelastic 
computations  in  3D,  in  contrast,  would  require  the  addition  of  roughly  30  memory 
variables  per  unit  cell  to  achieve  a  near-constant  Q  over  a  bandwidth  of  2  decades.  Thus, 
the  total  memory  requirement  would  be  roughly  quadrupled.  It  is  therefore  not  surprising 
that,  for  example,  none  of  the  3D  simulations  of  basin  response  to  earthquakes  done  to  date 
have  included  anelastic  memory  variables. 

We  present  a  method  that  approximates  constant  Q  to  very  high  accuracy,  over  2  to 
3  decades  of  frequency,  with  a  single  memory  variable  per  stress  component  per  unit  cell. 
The  method  is  most  efficient  in  higher  dimensions.  In  2D  the  method  achieves  a  4-fold 
reduction  in  memory  variable  storage  requirements  (to  achieve  Q  constant  over  a  given 
bandwidth),  compared  with  the  most  efficient  current  methods.  In  3D  the  method  provides 
an  8-fold  reduction.  The  operations  count  associated  with  the  memory  variable  evolution 
is  reduced  in  the  same  proportions. 

The  method  is  conceptually  simple  and  easily  implemented.  Given  an  N-pole 
description  of  the  anelastic  relaxation  function,  we  distribute  the  'N  memory  variables  over 
the  stress  node-points  of  the  computational  grid,  one  per  stress  component  per  unit  cell,  in 
such  a  way  that  spatial  averages  of  the  relaxation  function  remain  well  approximated.  This 
procedure  can  be  viewed  as  replacement  of  the  microscopic  relaxation  structure  of  the 
medium  with  a  coarse-grained  equivalent. 

The  following  section  reviews  the  operational  form  of  the  stress-strain  relationship 
of  linear  anelasticity,  its  representation  in  terms  of  a  relaxation  spectrum,  and  its 
approximation  via  first  order  differential  equations  governing  a  set  of  memory  variables. 
The  subsequent  section  analyzes  the  effect  of  coarse  graining  the  memory  variables.  The 
analysis  is  approximate,  and  is  intended  to  demonstrate  the  plausibility  of  the  technique  and 
to  indicate  how  the  evolutionary  equations  should  be  weighted  in  the  coarse-grained 
medium  to  achieve  the  desired  Q  value.  We  then  test  the  method  with  a  numerical 
application  to  constant  Q  simulation  in  3D.  Finally,  we  discuss  possible  extensions,  as 
well  as  limitations,  of  the  method. 

2.3  ANELASTICITY  AND  MEMORY  VARIABLES 

This  section  reviews  the  reduction  of  the  one-dimensional,  anelastic  stress-strain 
relationship  to  differential  form,  resulting  in  its  (approximate)  expression  in  terms  of  iV 
internal,  or  "memory",  variables.  In  higher  dimensions,  there  will  be  one  such  relationship 
for  each  independent  stress  component.  Adhering  to  the  notation  of  Day  and  Minster 
(1984),  we  write  the  one-dimensional  stress-strain  relation  as  a  convolution  integral  over 
the  step  response  M : 
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(1) 


c7(0  =  jM(r-f’)de(O, 

0 

where  a  is  the  stress  and  e  is  the  strain.  We  apply  the  Laplace  transform  in  s-multiplied 
form,  i.e., 


f{s)  =  s^f[t)e-‘*dt,  (2) 

0 

to  transform  the  stress-strain  relation  to  its  operational  form: 

a(s)  =  M(s)£(s).  (3) 

Note  that  the  operational  modulus  M  has  the  same  dimensions  as  the  step  response  M . 
The  unrelaxed  modulus  M^,  the  relaxed  modulus  Mr,  the  relaxation  of  the  modulus,  5M, 
and  the  normalized  relaxation  function  0,  are  given  by 


=  M(0)  =  M(oo)  , 

(4) 

Mg  =  M(oo)  =  M(0)  , 

(5) 

5M  =  M^-M,  , 

(6) 

Mit)  =  M„  +  5M(l)it)  . 

(7) 

We  represent  the  relaxation  function  in  terms  of  a  relaxation  spectrum  O, 


0{t)=  je-‘'^^{\nT)d(\nT)  ,  (8) 

resulting  in  the  following  integral  expression  for  the  operational  modulus: 


M  is)  =  -  5M  [  ^^^d(ln  t) 

L  ^1^+1  ,  (9) 

which  is  equivalent  to  Equation  8  of  Day  and  Minster  (1984).  We  use  (9)  and  the 
definition  of  (as  a  function  of  angular  frequency  co) 


Q-^(co)  = 


Im 

Miio)) 

Re 

M{ico) 

(10) 


to  get,  in  the  low  loss  approximation  (Q-^  «  1), 
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(11) 


,7  M\nr)d(lnt). 
cot;  + 1 


To  reformulate  (9)  in  differential  form,  we  approximate  the  integral  in  (9)  by  a 
numerical  quadrature,  i.e.,  a  sum  over  N  discrete  relaxation  times  r^. 


rOCnr) 
J  5T  +  1 


d(lnT)» 


A; 


tf^T,+l 


(12) 


where  the  A;  are  the  quadrature  weights.  Next,  we  substitute  the  approximation  into  (3) 
and  invert  the  Laplace  transform,  obtaining 


a(0  =  M„ 


£(0-S?.(0 


1=1 


(13) 


The  5;,  i=l,...N  are  then  internal,  or  memory,  variables  which  evolve,  respectively, 
according  to  the  N  first-order  differential  equations 


'■  dt 


+  C.(0  =  A;— e(0. 


(14) 


which  are  specified  by  the  2N  values  of  t,-  and  A, .  The  latter  are,  respectively,  the  poles 
and  residues  of  the  operational  modulus  (Day  and  Minster,  1984).  Each  stress  component, 
at  each  computational  node,  thus  has  associated  with  it  the  N  memory  variables  The 
resulting  low-loss  approximation  to  is  then 


(15) 


The  natural  criterion  on  which  to  judge  the  approximation  (12)  is  the  degree  to 
which  the  corresponding  approximation  (15)  to  Q'^(co)  matches  the  exact  value  (1 1).  Most 
seismic  applications  entail  approximating  a  Q'^{(o)  which  is  nearly  constant  over  a  broad 
frequency  range.  To  accomplish  this.  Day  and  Minster  (1984)  transformed  the  integration 
variable  in  (9)  to  -r'  and  applied  Gaussian  quadrature  to  derive  analytical  expressions  for 
the  poles  and  residues  T;  and  A,,  (and  showed  that  this  is  equivalent  to  a  Fade  approximant 
of  the  operational  modulus).  Subsequent  investigators  have  made  substantial 
improvements  by  using  numerical  optimization  to  determine  the  and  A,- ,  usually  leading 
to  roughly  uniform  distribution  of  the  poles  over  ln(T)  (e.g.,  Emmerich  and  Korn,  1987). 
For  practical  purposes,  about  1  pole  per  octave  is  sufficient  (Robertsson  et  al.,  1994; 
Blanch  et  al.,  1995)  to  approximate  a  frequency  independent  Q  function. 
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2.4  COARSE  GRAINING  OF  MEMORY  VARIABLES 
Outline  of  the  Method 

Our  objective  is  to  reduce  the  number  of  memory  variables  to  at  most  1  per  stress 
component  per  node,  without  sacrificing  accuracy  in  representing  Q(<o).  We  do  so  by 
introducing  a  single,  spatially  variable,  relaxation  time  t(x),  together  with  a  weighting 
function  w(x),  in  place  of  the  spectrum  of  relaxation  times  in  (8).  Then,  in  place  of  (13), 
we  write 


c7(x,r)  =  M„[e(x,0-C’U,0]. 

where  ^  satisfies 

T(x)^^^^  +  C(x,r)  =  w(x)^£(x,r).  (17) 

dt 

We  will  show  by  analysis  and  numerical  examples  that  this  method  can  be  made  highly 
successful  when  the  relaxation  time  function  t(x)  and  weight  function  w(x)  satisfy  two 
conditions.  First,  t(x)  and  w(x)  should  be  constructed  so  that  their  spatial  distributions 
approximate  the  target  relaxation  spectrum,  in  the  sense  described  in  the  next  paragraph. 
Second,  the  scale  length  of  their  spatial  variability  should  be  limited  to  a  few  times  the 
computational  cell  dimension,  as  shown  by  the  perturbation  analysis  in  the  next  section. 

Note  that  <l)(lnx)  in  Eqn.  8  is  non-negative  and  has  unit  area.  We  can  thus  think  of 
<E>(lnx)  as  a  probability  density  function,  in  which  case  the  integral  in  (9)  defining  the 
operational  modulus  would  be  the  expected  value  of  (5X+1)‘^  when  Inx  is  sampled  from 
the  distribution  d>.  It  is  natural  to  specify  x(x)  such  that  it  is  analogously  distributed,  in  the 
following  sense.  Choose  x(x)  to  be  a  measurable  function  on  domain  X,  and  define  the  set 
A  by 

A  =  {x  :  Inx(x)  <a}  (18) 

where  a  is  a  real  number.  We  introduce  a  measure  v,  such  that  the  measure  of  a  subset  E  of 
X  is  defined  by 

v(£)  =  |w(x)<iV,  (19) 

E 

where  h'(x)  is  a  Lebesgue-measurable,  non-negative  function  having  unit  mean  value  on 
domain  X.  This  means  that 

v(X)  =  m(X),  (20) 
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where  m  is  the  Lebesgue  measure  of  X  (i.e.,  the  volume  of  the  domain  of  interest).  Then 
we  would  specify  x(x)  such  that  the  measure  of  A,  v(i4),  satisfies 


J<D(lnT)rf(lnT), 


(21) 


for  all  real  a.  In  other  words,  given  a  w(x)  satisfying  (20),  we  specify  t(x)  such  that,  if 
spatial  points  were  sampled  randomly,  with  probability  weighted  by  w(x),  the  resulting  Inx 
samples  would  have  distribution  O.  In  practice,  of  course,  we  will  only  need  to  specify  w 
and  X  at  a  discrete  set  of  computational  nodes. 

There  are  actually  two  distinct  scales  of  spatial  variation  implicit  in  Equation  17. 
The  micro-scale  (sub- wavelength)  variation  which  we  have  just  discussed,  and  which  is 
designed  to  simulate  the  relaxation  spectrum,  is  indicated  explicitly  by  the  notation  in  (17). 
Equation  17,  however,  does  not  explicitly  note  the  spatial  variation  of  the  factor 
though  it  too  may  vary  spatially.  This  factor,  however,  scales  the  macroscopic  Q,  and  its 
variation  will  usually  have  a  scale  length  larger  than  that  of  x  and  w. 

We  will  refer  to  the  conventional  description  of  anelasticity,  as  given  by  Eqn  9,  as 
the  "relaxation-spectrum"  representation.  We  will  call  the  method  based  on  a 
heterogeneous  relaxation  time,  as  given  by  Eqns  16  and  17,  the  "coarse-grain" 
representation,  in  accordance  with  the  physical  analogy  described  below. 


Physical  analogue 

Some  motivation  for  this  approach  and  terminology  may  be  found  in  the  physical 
origin  of  the  relaxation  spectra  of  earth  materials.  The  characteristic  attenuative  behavior  of 
rock  results  from  the  fact  that  rock  is  a  heterogeneous  aggregate.  Relaxation  mechanisms 
corresponding  to  a  broad  range  of  scale  lengths  and  activation  energies  are  present  on  the 
microscopic  level  (see,  e.g..  Minster,  1980).  When  the  material  is  deformed  over 
wavelengths  large  compared  to  the  scale  of  the  heterogeneities,  these  relaxation 
mechanisms  are  activated  concurrently,  giving  rise  to  a  broad  spectrum  of  relaxation  times. 
Conventionally,  these  relaxations  are  subsumed  into  a  constitutive  equation  describing  the 
mean  behavior  of  the  aggregate.  Our  method,  in  a  sense,  reverses  this  synthesis.  We 
revert  to  the  more  fundamental,  heterogeneous  description  of  the  relaxations,  albeit  in  an 
artificial,  coarser-grained  aggregate,  with  "grain  size"  equal  to  the  computational  unit  cell. 
In  this  coarse  aggregate,  each  computational  unit  cell  has  only  a  single  relaxation  time. 

In  discrete  time-stepped  numerical  wave  propagation  computations  (e.g.,  finite 
difference,  finite  element,  pseudospectral),  the  minimum  wavelength  which  can  be 
computed  accurately  necessarily  exceeds  several  times  the  cell  size  of  the  grid.  Thus,  it  is 
plausible  that  we  can  coarse-grain  the  medium  in  such  a  way  that  its  attenuative  properties 
will  be  preserved  for  all  wavelengths  of  computational  interest.  The  following  analysis  of 
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the  coarse-grain  representation  by  means  of  perturbation  theory  suggests  that  this  is  indeed 
the  case. 

Analysis  of  the  method 

In  the  coarse-grain  representation  of  the  medium,  the  effective  attenuation  can  no 
longer  be  regarded  as  a  purely  local  process  with  properties  derivable  from  constitutive 
equations  alone.  To  understand  the  nature  of  the  approximation,  it  is  necessary  to 
introduce  the  elastodynamic  equations  into  the  analysis.  We  begin  with  the  well-known 
analysis,  via  perturbation  theory,  of  anelastic  attenuation  in  terms  of  an  eigenfunction 
expansion  of  the  elastodynamic  field.  We  simplify  the  analysis,  without  sacrificing  any 
useful  insights,  by  restricting  consideration  to  wavefields  in  a  bounded  volume  (and  by 
neglecting  issues  raised  by  eigenfrequency  degeneracy). 

In  a  finite  volume  V  with  homogeneous  boundary  conditions,  all  solutions  to  the 
elastodynamic  equations  (without  anelastic  losses)  can  be  written  as  linear  combinations  of 
u„(x)e"””‘,  /i=l,2  ...,  where  is  the  (normalized)  eigenfunction  belonging  to  the 
eigenvalue  co^  of  the  operator  given  by 

H?(u)  =  p->(C^U„),j.  (22) 

In  (22),  u(x)c‘“‘  is  the  displacement  field,  p  (x)  the  density,  and  C®(x)  the  elastic  tensor. 
The  eigenvalues  are  positive  and  the  eigenvectors  are  orthogonal  under  the  inner  product 
defined  by 


(f,g)  =  Jf(x)-g*(x)p(x)dy.  (23) 

V 

The  introduction  of  weak,  anelastic  attenuation  has  the  effect  of  adding  a  small 
additional,  frequency  dependent  term  H'  to  H°.  The  new  operator  H  is  then 

H  =  H°-hH',  (24) 


where 


h;(u)=p-'(c;,u„).;  .  (25) 

Each  component  of  C'  has  the  form  of  the  second  term  on  the  right  hand  side  of  (9), 
evaluated  at  5=/co .  The  resulting  perturbation  to  the  eigenvalue  col  can  be  estimated  from 
perturbation  theory.  The  first-order  correction  5col  to  the  eigenvalue  is 

5a)/  =  (H'(u„),u„),  (26) 
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where  Un  is  the  normalized,  unperturbed  eigenfunction  of  The  temporal  Q'^  (Aki  and 
Richards,  1980)  of  the  mode,  Q~' ,  can  be  expressed  in  terms  of  to  first  order  in  the 
components  of  C',  as 


(2;>=2a);4m(5a)J, 

=  (i);2Im(H'(uJ.u«) 


(27) 


To  simplify  the  subsequent  discussion,  we  now  specialize  to  the  case  of  an 
isotropic  medium,  though  this  assumption  is  not  essential  to  the  analysis.  To  begin  with, 
we  consider  the  relaxation-spectrum  representation  of  the  medium.  We  take  a  normalized 
relaxation  spectrum  O  which  is  independent  of  position  x  (though  this  assumption  is  not 
essential  either).  Then  we  can  find  H'  from  (9),  replacing  M  and  6M  by  their  bulk  and 
shear  modulus  equivalents: 

H,(u)  =  -p-‘(5.s„«,  +  2«K)J^dlnr.  (28) 

where  e  is  the  infinitesimal  strain  tensor  and  e'  its  deviatoric  part,  and  k  and  |X  denote, 
respectively,  bulk  and  shear  moduli.  Then,  substituting  (28)  in  (27),  and  applying  the 
divergence  theorem  and  the  boundary  conditions,  we  obtain 

~  r 

Qn'  =  J ^<0.  (t^)^(ln  t)d In T 

-«o  ^ 

V 


^^uKn)ik^(n)}k 


W,  (29) 


where  £nand  e'n  are  strain  tensors  associated  with  the  eigenvector  u„,  and  i7o>('C)  is 
defined  by 


X(0 

1  +  T^Q}^  ' 


(30) 


Recalling  our  reinterpretation  of  <I>(lnx)  as  a  probability  density  function,  we  can 
interpret  the  first  integral  in  (29)  as  the  expected  value  of  (t).  Denoting  this  expected 
value  by  and  also  introducing  and  £nn  to  denote,  respectively,  the  bulk  and 

shear  strain  energy  densities  of  mode  n,  we  can  write  (29)  in  the  form 


Now  we  are  ready  to  introduce  the  coarse-grain  representation,  as  given  by  Eqn  16- 
21.  Therefore,  we  drop  the  integral  over  the  relaxation  spectrum,  i.e.,  our 

representation  of  H',  and  instead  allow  the  relaxation  time  x  to  be  a  function  of  position. 
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x(x),  as  described  earlier.  The  result  is  that  in  (31)  is  replaced  by  w(\)D^Jr(x)). 
We  denote  by  Q,  the  resulting  approximation  to  the  temporal  Q  of  mode  n.  Then,  Q*'  i^ 


EA 


Hx)dV. 


(32) 


We  can  assess  the  coarse-grain  representation  by  evaluating  the  accuracy  of  (32)  as 
an  approximation  to  (31),  since  the  latter  is  the  relaxation- spectrum  representation  of  the 
eigenfunction  attenuation.  The  comparison  is  more  informative  with  the  further 
simplification  that  5k/Ku  and  5|i/pu  are  spatially  uniform.  From  its  method  of  construction, 
we  know  that  the  volume  average  of  w(x)D^^  ('^(*))  is  also  equal  to  the  expected  value 
(d^  ^(see  Eqns  18-21).  We  write  and  for  the  corresponding  volume  averages  of 
EnK  and  E^i,  and  (31)  simplifies  to 


er=2a)rv 


-2l 

n 


ir  ^  II 

\^u  hu 


KY 


(33) 


The  corresponding  expression  for  attenuation  of  the  eigenfunction  in  the  coarse-grain 
representation  is 


q:'=2(o: 


— j£„(x)D.,  +  2^  J  £„(x)D.,  (r{x)Mx)<<V' 


=  2.-' +  2®.-’ 


r 


Sk 


5n 


— C,  +  2^C„ 


(34) 


The  Ck  and  in  (34)  are  just  the  spatial  correlations  of  w(x)D^^  ('r(x))  with  the  bulk  and 
shear  strain  energy  densities,  respectively: 

c«  =|[£.«(x)-£.„lH'(x)D._(T(x))-:jD..]dl',  (35) 


where  M  refers  to  either  k  or  |i,  and  wDm,  is  the  volume  average  of  w{\)D^^  ('r(x))-  We 
have  used  the  equality  between  wDm,  and  ^  in  deriving  (34)  and  (35).  Eqn.  35  states 
that  the  error  introduced  by  the  coarse-grain  representation,  relative  to  the  relaxation- 
spectrum  representation,  is  proportional  to  the  degree  of  spatial  correlation  between  the 
modal  strain  energy  density  and  the  function  w(\)D^^  (■r(x)). 

We  should  be  able  to  minimize  the  correlation  of  w(\)D^^  i'^W)  with  all  modes  of 
practical  interest  by  making  x(x)  and  w(x)  periodic  in  space,  with  wavelength  short 
compared  with  the  minimum  wavelength  at  which  the  numerical  wave  propagation 
computation  (e.g.,  finite  difference,  etc.)  is  considered  accurate.  We  assess  this  strategy 
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by  examining  (34)  and  (35)  in  the  simplified  case  of  plane  wave  excitation  of  a  uniform,  ID 
acoustic  medium. 

In  this  special  case,  the  eigenfunction  has  the  form 

ujx)  =  5cos(<u,Jc/a) ,  (36) 

(apart  from  an  irrelevant  phase  shift),  where  a  is  the  wavespeed  and  is  a  constant.  The 
deviation  of  the  strain  energy  density  from  its  mean  value  of  af  is 

(x)  -E^=~k^  (co^bI  a)^  cos(2tu,x/a) ,  (37) 

which  has  1/2  the  wavelength  of  the  eigenfunction.  We  will  assume  that  »v(x)  and  x(x)  are 
periodic,  with  period  L,  i.e.,  t(x)  =  x{x-mL)  for  integer  m,  and  similarly  for  w.  Then  the 
function  w(\)D^  (t(x))  is  also  periodic  with  period  L,  and  therefore  has  a  Fourier  series 
representation,  the  zeroth  order  term  of  which  is  just  wDa. ,  so  that  the  deviation  of 
('^(x))  front  its  mean  is 

oo  OS 

M/(x)D^^(T(x))-M^a).  ='^ajCOsi27j:jxlL)  +  '^bjSmi27gxfL).  (38) 

y=i 

Combining  (37)  and  (38)  with  (35),  we  find  that  the  correlation  Ck  is  only  appreciable 
when  the  condition 


^27ca^ 


V  y 


=  2L, 


(39) 


is  satisfied,  where  27ca/C0n  is  the  wavelength.  Therefore,  the  largest  value  of  wavelength 
for  which  the  error  in  the  coarse-grain  representation  is  significant  is  2L .  In  other  words, 
the  estimated  error  is  only  significant  when  the  wavelength  is  less  than  or  equal  to  twice  the 
periodicity  of  the  relaxation  time  function  t(x).  Note,  by  the  way,  that  (39)  is  just  the 
Bragg  condition  for  backscatter  at  normal  incidence  to  the  scattering  plane. 


Implementation 

Eqn.  39  implies  a  bandwidth  limitation  in  the  coarse-grain  method.  Our  strategy  for 
its  implementation  is  to  match  this  short-wavelength  cutoff  to  the  short-wavelength  limit 
already  present  in  discrete  numerical  methods.  We  make  w(x)  and  t(x)  periodic,  with 
period  equal  to  twice  the  computation  unit  cell  dimension.  According  to  (39),  this  period- 
two  coarse-graining  will  not  result  in  significant  error  for  wavelengths  longer  than  4  unit¬ 
cell  dimensions.  For  relatively  low-order  finite  difference  methods,  including  the  widely 
used  fourth  order  staggered  grid  methods  (e.g.,  Levander,  1988),  wavelengths  shorter  than 
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4  unit-cell  dimensions  suffer  strong  numerical  dispersion  and  are  therefore  already  subject 
to  large  errors.  Thus,  any  spurious  effects  of  period-two  coarse-graining  should  be 
isolated  at  wavelengths  which  are  already  outside  the  usable  bandwidth  of  the  wave 
propagation  calculations.  High-order  finite  difference  methods  and  pseudospectral 
methods  are  less  restricted  by  numerical  dispersion,  but  suffer  other  limitations  associated 
with  the  representation  of  material  interfaces  (Witte  and  Richards,  1990)  when  the 
wavelength  is  shorter  than  4  times  the  cell  dimension,  and  the  period-two  coarse-grain 
method  may  turn  out  to  be  satisfactory  in  those  applications  as  well. 

In  3D,  this  2-unit-cell  periodicity  allows  for  8  different  values  of  t  at  grid  nodes. 
To  be  specific,  if  this  chosen  set  of  discrete  relaxation  times  is  and  unit  cells 

are  indexed  by  integers  p,q,r,  such  that  stress  nodes  are  centered  at  x  =  p&i+qai  -i-raj  , 
then 

T(x)  =  Tj^.^^oj2+2(?mod2)+4(rmod2)» 

and  wCx)  has  a  similar  discrete  representation 


w(x) 

^l+pmod2+2(<ymod2)+4(rmod  2)  * 


Since  this  approach  only  permits  a  discrete  set  of  relaxation  times  (8  in  3D)  to  be 
realized,  implementation  of  criterion  (21)  requires  first  developing  a  discrete  approximation 
to  the  continuous  relaxation  spectrum,  i.e., 

N 

C>(lnT)  =  ^A,5(lnT-lnT,),  (42) 

i=l 


as  in  (12).  This  is  the  same  step  required  in  the  generation  of  a  conventional  memory 
variable  representation,  and  any  of  the  optimization  methods  proposed  in  the  literature 
(e.g.,  Emmerich  and  Korn,  1987;  Carcione  et  al.,  1988;  Witte  and  Richards,  1990;  Blanch 
et  al.  1995)  for  generating  conventional  memory  variable  representations  can  be  used  at  this 
point.  Once  the  continuous  spectrum  is  replaced  by  the  discrete  spectrum,  the  conditions 
defined  by  (18-21)  will  be  satisfied  if  the  Wi  are  identified  with  the  quadrature  weights  Xi, 
but  normalized  to  have  unit  volumetric  mean.  Since  the  normalization  of  O  implies  that  the 
Xi  sum  to  1,  we  then  have 

Wi  =  NXi,  (43) 


where  is  8  in  3D  or  4  in  2D.  Fig.  1  illustrates  this  period-two  scheme. 

Eight  relaxation  times  should  be  sufficient  to  approximate  a  frequency-independent 
QioS)  with  very  high  accuracy  over  roughly  3  decades  in  frequency.  In  2D,  the  2-unit-cell 
periodicity  allows  for  4  distinct  values  of  x,  adequate  to  approximate  frequency- 
independent  Qia)  over  one  and  a  half  decades  or  more.  Compared  with  a  conventional 
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memory  variable  method,  i.e.,  one  in  which  the  memory  variables  are  not  coarse  grained, 
the  coarse-grain  method  with  2-unit-cell  periodicity  can  approximate  Q(Q))  with  the  same 
accuracy,  but  reduces  the  memory  requirement  by  a  factor  of  8  in  3D  and  a  factor  of  4  in 
2D. 


2.5  NUMERICAL  EXAMPLE 

The  above  analysis  is  approximate,  and  a  numerical  test  is  necessary  to  establish  the 
practical  accuracy  of  the  coarse-grain  method.  In  this  section,  we  test  its  period-two 
implementation  for  the  case  of  a  plane  wave  in  a  uniform,  3D  acoustic  medium,  with  Q 
approximately  frequency-independent.  We  solve  the  3D  acoustic  equations  using  a 
staggered  grid  method,  fourth  order  in  space  and  second  order  in  time  (Graves,  1996, 
describes  the  corresponding  algorithm  for  the  elastodynamic  equations).  Eqn  17  is  time 
stepped  using  the  formula 

C(x,r  +  &)  =  +  w(x)— i[£(x,r)  -i-  e(x,t+ &)|l  -  (44) 

1C,. 


A  relaxation  spectrum,  constant  on  its  support  interval  (InXmdn'^M)  yields  a  good 
approximation  to  a  frequency  independent  Q  for  TM'^«co«Xm‘^  (this  is  a  change  of 
notation  from  Day  and  Minster,  1984,  in  that  Xm  and  Xm  correspond  to  their  Xi  and  X2).  If 
©0  is  a  reference  frequency  near  the  center  of  the  absorption  band,  then  =  Q~^(cOo)  is 
given  approximately  by 


nSic 


2  k.. 


In 

fx  ^ 

-K 

8k,  1  \ 

+— In  ©o^«) 

J 

K  \ 

(45) 


To  approximate  this  absorption  band,  we  take  h'(x)  equal  to  1  everywhere,  and  prescribe 
x(x)  according  to  (40),  with  the  InXk  uniformly  spaced  over  the  interval  (lnXm,lnXM): 

In  T*  =  In  T„  +  ^(In  -  In  rj.  (46) 

The  plane  wave  is  given  an  initial  displacement  pulse  shape  uo(t),  where 

which  is  a  broadband,  minimum  phase  pulse  with  its  spectral  comer  at  frequency  1/T.  In 
the  cases  shown  below,  T  is  set  equal  to  lOtt  times  the  propagation  time  across  the  unit 
cell,  so  that  the  comer  frequency  corresponds  to  a  wavelength  of  about  30  grid-cells.  The 
propagation  direction  is  parallel  to  one  of  the  lattice  vectors  of  the  grid. 
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We  have  set  the  ratio  XM/Xm  to  KH,  with  the  objective  of  producing  a  nearly 
frequency  independent  Q  for  IOxm'^  <  to  <  O.lXm’^  Co  is  defined  at  reference  frequency 
coq,  which  we  set  equal  to  the  geometrical  mean  of  xm*^  and  Xm'^  We  perform  calculations 
for  Co  of  100  and  20.  The  lower  cutoff  Xm  is  set  equal  to  (27t)'i  times  the  propagation  time 
across  the  unit  cell. 

Fig.  2  shows  the  pulse  for  the  case  Co  equal  to  100,  after  it  has  propagated  a 
distance  of  26  times  the  comer  wavelength.  Fig.  2a  shows  the  entire  pulse,  with  the  later, 
low-amplitude  parts  of  the  waveform  magnified  by  factors  of  100  and  1000.  Fig.  2b 
shows  the  high  amplitude  part  of  the  pulse  on  an  expanded  time  scale.  For  comparison, 
both  parts  of  Fig.  2  also  show  the  exact  solution,  (based  on  the  frequency-independent  C 
model  of  Kjartansson,  1979).  It  is  evident  that  the  coarse-grain  method  reproduces  the 
analytical  solution  with  very  high  accuracy.  The  small,  negative  precursor  in  the  numerical 
solution  is  a  consequence  of  ordinary  numerical  dispersion  associated  with  the  fourth  order 
finite  difference  method,  and  is  unrelated  to  the  anelastic  relaxation  model. 

Fig.  3  shows  the  corresponding  results  for  Co  equal  to  20.  As  in  the  previous  case, 
the  numerical  and  exact  solutions  are  compared  in  the  figure,  with  the  late  portion 
magnified  in  Fig.  3a  and  the  early  portion  expanded  in  Fig.  3b.  The  numerical  solution 
reproduces  both  the  attenuation  and  physical  dispersion  of  the  exact  solution  with  very  high 
accuracy.  Fig.  4  compares  the  analytical  and  numerical  solutions  again,  this  time  after 
convolution  with  a  Ricker  wavelet  with  dominant  frequency /r  equal  to  twice  the  comer 
frequency.  Thus,  the  pulse  has  propagated  52  times  its  dominant  wavelength.  Again, 
numerical  and  analytical  solutions  are  virtually  indistinguishable,  even  at  late  times. 

We  can  determine  an  apparent  C.  denoted  by  Qico),  from  the  numerical  solutions. 
We  do  this  using  spectral  ratios  of  the  pulse  at  successive  locations,  i.e., 

Q-U(0)  =  + 

coaAx 

where  Ax  is  the  propagation  distance  between  locations  and  a  is  the  wavespeed.  This 
operational  definition  of  Q  differs  from  that  of  Eqn  10,  but  at  large  values  of  Q,  the 
difference  is  negligible  (see,  e.g.,  O'Connell  and  Budiansky,  1978).  The  difference  can 
become  significant  for  small  Q,  and  this  difference  will  become  evident  in  our  case  when 
j2o  equals  20. 

Fig.  5  shows  Qico)  for  the  numerical  solution  in  the  case  Qo  =  100.  Four  estimates 
are  shown,  corresponding  to  spectral  ratios  at  successive  200  grid-cell  intervals.  The 
Qico)  curves  are  indistinguishable  from  one  another  over  the  entire  range  of  frequencies 
from  1(P  grids  per  wavelength  to  4  grids  per  wavelength.  Thus,  the  coarse-grain  method 
yields  a  spatially  homogeneous  attenuative  behavior  over  the  full  range  of  usable 
frequencies  present  in  the  finite  difference  computations.  The  variability  of  the  curves  at 
wavelengths  longer  than  10^  grid-cells  is  simply  a  result  of  the  shorter  pulse  durations 
available  for  Fourier  analysis  at  greater  distances  (due  to  later  arrival  times),  10^  grids  per 
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wavelength  corresponding  to  a  period  equal  to  the  total  post-arrival  duration  of  the 
computation  at  the  most  distant  location.  Likewise,  4  grids  per  wavelength  is  well  below 
the  practical  limit  imposed  by  numerical  dispersion  in  the  fourth  order  method  used  here 
(Levander,  1988). 

In  addition  to  indicating  a  spatially  homogeneous  behavior.  Fig.  5  shows  that  the 

^  _ 

method  is  highly  successful  in  producing  a  frequency-independent  Q(ci)).  The  dotted  box 
encloses  the  usable  computational  bandwidth  and  indicates  the  6%  tolerance  interval  about  a 
Q  of  100.  For  the  entire  computational  bandwidth,  Q  lies  within  6%  of  the  target  value  of 
100.  Fig.  6  shows  the  corresponding  results  for  (2o  =  20.  In  this  case,  the  difference  in 
definition  between  Q  and  Q  is  significant,  and  when  Q  (as  defined  by  Eqn  10)  is 
frequency  independent  and  equal  to  20,  Q  (as  defined  by  Eqn  48)  should  be 
approximately  19,  as  can  be  verified  analytically  using  expressions  in  Kjartansson  (1979). 
As  Fig.  6  indicates,  this  target  Q  is  again  achieved  within  6%  tolerance  over  virtually  the 
entire  usable  computational  frequency  band.  For  the  2  decades  from  5  grids  per 
wavelength  to  500  grids  per  wavelength,  the  variation  in  Q  is  less  than  ±3%  in  both  the 
j2o  =  20  and  Qo  =  100  cases. 

At  wavelengths  shorter  than  4  grid  points,  as  predicted  by  the  perturbation  analysis, 
backscatter  due  to  coarse-graining  adds  to  the  numerical  errors  which  are  already  present  in 
the  form  of  numerical  dispersion.  For  a  low-order  finite  difference  method  (such  as  the 
fourth  order  method  tested  here),  wavelengths  shorter  than  4  grid  points  are  usually 
considered  to  lie  outside  the  usable  computational  bandwidth,  as  a  result  of  ordinary 
numerical  dispersion  alone.  Whenever  that  is  the  case,  the  anelastic  coarse-graining 
imposes  no  additional  limitation  on  usable  bandwidth.  On  the  other  hand,  higher-order 
finite  difference  and  pseudospectral  methods  may  have  sufficiently  accuracy  for  some 
applications  at  wavelengths  shorter  than  4  grid  points.  When  that  is  the  case,  the  coarse- 
grain  method  will  not  be  appropriate. 

2.6  DISCUSSION 

In  the  numerical  example,  the  period-two  coarse-grain  method  achieved  a 
frequency-independent  Q,  to  high  accuracy,  over  more  than  2  decades  in  frequency,  using 
a  single  memory  variable  per  stress  component  per  node  .  This  accuracy  was  achieved 
without  any  attempts  to  optimize  the  frequency  distribution  or  weighting  of  the  relaxation 
times  in  (46).  In  fact,  we  did  not  exploit  the  flexibility  provided  by  the  weight  function  w 
at  all,  simply  setting  it  to  1  in  the  example. 

The  weight  function  can  be  exploited  in  several  different  ways  to  optimize  and/or 
generalize  the  coarse-grain  method.  For  example,  if  further  accuracy  were  desired  in  the 
approximation  of  a  frequency-independent  Q,  any  of  the  optimization  methods  which  have 
been  proposed  in  the  literature  for  generating  conventional  memory  variable  representations 
could  be  used  in  conjunction  with  period-two  coarse  graining  (e.g.,  Emmerich  and  Korn, 
1987;  Carcione  et  al.,  1988;  Witte  and  Richards,  1990;  Blanch  et  al.  1995).  The  first  step 
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in  doing  so  would  be  the  conventional  one:  choose  relaxation  times  Xi  and  quadrature 
weights  Xi  (8  of  each  for  3D,  4  of  each  for  2D)  to  optimize  the  approximation  of  (11)  by 
(15).  Then  those  relaxation  times  would  be  identified  with  the  discrete  form  of  x(x),  Eqn 
40;  the  quadrature  weights,  after  being  normalized  to  unit  mean,  would  be  identified  with 
the  discrete  form  of  w(x),  Eqn  41. 

A  second  application  of  the  weight  function  would  be  to  approximate  some  other 
desired  frequency  dependence  of  Q,  a  power  law,  for  example.  As  in  the  case  of  frequency 
independent  Q,  any  of  the  memory  variable  optimization  techniques  in  the  literature  could 
be  used  to  find  the  8  relaxation  times  and  quadrature  weights,  which  are  then  identified 
with  the  coarse-grain  relaxation  times  and  (after  normalization)  weights. 

A  third  application  of  the  weight  function  would  be  to  further  reduce  memory 
requirements  in  instances  where  the  accuracy  possible  with  8  relaxation  times  is  redundant. 
In  some  seismic  exploration  problems,  for  example,  a  fairly  narrow-band  representation  of 
Q  may  be  acceptable  in  practice.  Blanch  et  al.  (1995)  have  argued  that  2  properly  chosen 
relaxation  times  provide  sufficient  constant-^  bandwidth  (Q  constant  within  8%  over 
roughly  one  decade  in  frequency)  for  many  practical  purposes,  especially  given  the  limited 
problem  sizes  which  are  practical  with  current  computers  (Robertsson  et  al.  (1994)  have 
even  found  a  single  relaxation  time  to  provide  an  adequate  approximation  in  some  narrow- 
band  applications  in  high  resolution  reflection  seismology).  In  cases  where  the  two- 
relaxation-time  approximation  is  acceptable,  for  example,  we  could  assign  these  2  values  to 
Ti  and  X2  in  (40),  assign  values  to  wi  and  wz  in  (41)  such  that  wi  -h  w2  equals  8  (in  3D;  or 
4  in  2D),  and  set  the  remaining  wj  to  zero.  This  strategy  would  eliminate  the  need  to  store 
any  memory  variables  at  all  at  75%  of  the  nodes  in  the  grid,  and  only  1  per  node  at  the 
remaining  25%,  thereby  reducing  the  memory  variable  storage  by  a  factor  8  (or  a  factor  of 
4  in  2D),  relative  to  a  conventional  2  memory-variable  method.  Obviously,  then,  with  the 
aid  of  the  weight  function,  any  conventional  memory  variable  scheme  could  be  coarse 
grained  to  achieve  an  8-fold  memory  reduction  in  3D,  or  4-fold  reduction  in  2D. 

The  above  numerical  example  was  restricted  to  the  acoustic  wave  equation.  There 
is  no  apparent  reason,  however,  that  the  method  should  not  work  equally  well  for  the 
elastodynamic  case.  In  fact,  it  is  for  the  later  case  that  the  method  is  likely  to  be  most 
useful,  because  the  presence  of  6  stress  components  in  3D  elastodynamics  places  very  large 
memory  requirements  on  the  conventional  memory  variable  method.  Furthermore,  the 
advantages  of  the  coarse-grain  method  will  become  even  more  important  as  higher 
computing  speeds  become  available.  The  speed  increases  will  permit  3D  elastodynamic 
simulations  to  be  computed  with  ever  greater  usable  bandwidth,  and,  as  computational 
bandwidth  increases,  memory  variables  constitute  a  progressively  larger  proportion  of  the 
total  computational  memory  requirement. 

We  have  established  the  viability  of  the  coarse-grain  method  for  the  fourth-order 
staggered  grid  finite  difference  method,  and  it  should  work  equally  well  with  other  low- 
order  finite  difference  or  finite  element  methods.  It  has  yet  to  be  established,  however. 
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whether  the  coarse-grain  method  will  work  effectively  with  higher-order  finite  difference 
methods  or  with  pseudospectral  methods.  The  main  issue  is  likely  to  be  whether  coarse- 
graining  entails  any  loss  of  usable  bandwidth  when  applied  with  those  numerical  methods. 
The  staggered-grid  pseudospectral  methods,  in  particular,  may  be  questionable  candidates 
for  the  coarse-grain  method.  In  pseudospectral  methods,  the  spatial  differentiation 
operators  are  exact  if  the  medium  is  homogeneous  (and  the  wavefield  spatially  periodic  and 
sufficiently  band-limited).  In  that  case,  the  only  source  of  numerical  dispersion  arises  firom 
approximations  associated  with  the  time  integration.  In  the  case  of  the  staggered  grid 
pseudospectral  methods,  numerical  dispersion  can  be  kept  very  small  even  at  the  spatial 
Nyquist  wavenumber,  i.e.,  down  to  only  2  nodes  per  wavelength.  On  the  other  hand, 
when  discrete  interfaces  or  sharp  transitions  in  wavespeed  are  present  in  the  model,  3  to  4 
nodes  per  wavelength  are  required  in  order  to  accurately  model  the  reflection  and 
transmission  behavior  (Witte  and  Richards,  1990).  In  cases  where  the  practical  limit  of  the 
pseudospectral  method  is  determined  to  be  4  nodes  per  wavelength  or  greater,  the  coarse- 
grain  method  should  be  effective,  though  this  conjecture  will  have  to  be  tested  numerically. 
If  so,  the  method  will  provide  the  same  degree  of  memory  savings  that  it  provides  in  the 
case  of  low-order  finite  difference  methods. 

2.7  SUMMARY 

We  have  developed  a  new  method  for  introducing  anelastic  attenuation  into  time- 
stepped  numerical  methods  for  solving  the  equations  of  acoustics  and  elastodynamics.  The 
method  is  designed  to  reduce  drastically  the  memory  requirements  of  anelastic  modeling.  It 
takes  as  its  starting  point  an  7V-pole  (i.e.,  A^-relaxation  time)  approximation  to  the  relaxation 
spectrum.  Then,  the  N  corresponding  memory  variables  are  distributed  over  the  stress 
node-points  of  the  computational  grid,  one  per  stress  component  per  unit  cell.  An  analysis 
via  perturbation  theory  shows  that  this  single  memory-variable  method  will  reproduce  the 
attenuative  behavior  of  the  N  memory-variable  model  provided  i)  the  weights  assigned  each 
of  the  N  relaxations  (corresponding  to  the  residues  in  the  TV-pole  representation)  are 
rescaled  to  have  unit  volumetric  mean,  and  provided  the  relaxation  times  are  distributed  in  a 
periodic  array,  with  period  less  than  half  the  minimum  wavelength  to  be  modeled.  We 
make  N  equal  to  8  for  3D  problems  and  4  for  2D  problems,  and  distribute  the  memory 
variables  so  that,  for  each  stress  component,  the  relaxation  time  array  has  a  spatial  period, 
in  each  direction,  equal  to  twice  the  unit-cell  dimension.  This  period-two  coarse-grain 
method  is  accurate  for  wavelengths  exceeding  4  unit-cell  dimensions. 

Accuracy  of  the  method  has  been  verified  in  numerical  examples  for  the  3D  acoustic 
wave  equation.  The  test  problems  were  solved  using  a  fourth  order  staggered  grid  finite 
difference  method.  In  the  numerical  examples,  a  frequency-independent  Q  was  achieved 
within  a  3%  tolerance  over  the  2  decades  in  frequency  from  5  grids  per  wavelength  to  500 
grids  per  wavelength,  using  only  a  single  memory  variable  per  stress  component  per  node. 
Straightforward  generalizations  of  the  method  can  be  used  to  approximate  specific 
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frequency-dependent  Q  models  such  as  power  laws,  or  to  further  reduce  memory 
requirements. 

We  expect  the  method  to  work  equally  well  for  low-order  finite  difference  or  finite- 
element  solutions  in  elastodynamics.  For  higher-order  finite  difference  methods,  as  well  as 
for  pseudospectral  methods,  the  utility  of  the  coarse-graining  method  will  depend  upon  the 
limiting  wavelength  at  which  the  underlying  numerical  method  has  acceptable  accuracy.  In 
those  numerical  methods  in  which  the  usable  bandwidth  is  restricted  to  wavelengths  longer 
than  4  unit-cell  dimensions,  the  coarse-grain  method  should  be  accurate  and  efficient. 

In  those  cases  where  it  is  appropriate,  the  method  reduces  storage  requirements  for 
memory  variables  by  a  factor  of  8  for  3D  problems  and  by  a  factor  of  4  for  2D  problems. 
As  computer  speeds  increase,  larger  problems,  and  therefore  higher  computational 
bandwidths,  will  become  feasible.  As  computational  bandwidth  increases,  memory 
variables  will  account  for  an  increasing  proportion  of  total  storage  requirements,  putting  a 
further  premium  on  the  concision  achieved  by  the  coarse-grain  method. 
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Figure  2  (a).  Broadband  pulse,  for  the  case  of  Qq  equal  to  100,  after  propagating  a  dist^ce 
of  26  times  its  comer  wavelength.  The  entire  pulse  is  shown,  with  low-amplitude  portions 
magnified  by  factors  of  100  and  1000.  Time  has  been  scaled  by  the  source  comer 
frequency  1/T. 
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Figure  2  (b).  Broadband  pulse,  for  the  case  of  Qq  equal  to  100,  after  propagating  a 
distance  of  26  times  its  corner  wavelength.  The  initial  portion  is  shown  on  an 
expanded  time  scale.  Time  has  been  scaled  by  the  source  comer  frequency  l/T. 


Normalized  Displacement 


Figure  3  (a).  Broadband  pulse,  for  the  case  of  Qo  equal  to  20,  after  propagating  a  distance 
of  26  times  its  comer  wavelength.  The  entire  pulse  is  shown,  with  low-amplitude  portions 
magnified  by  factors  of  10  and  100. 
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Figure  3  (b).  Broadband  pulse,  for  the  case  of  Qq  equd  to  20,  after  propagating  a  distance 
of  26  times  its  comer  wavelength.  The  initial  portion  is  shown  on  an  expanded  time  scale. 


Figure  4  (a).  Narrow-band  pulse,  for  the  case  of  Qo  equal  to  20,  after  it  has  propagated  a 
distance  of  52  times  its  dominant  frequency.  The  pulse  was  constructed  by  convolving  the 
broadband  pulse  with  a  Ricker  wavelet.  The  latter  has  a  dominant  frequency^  equal  to 
twice  the  comer  frequency  of  the  initial  broadband  pulse.  The  entire  pulse  is  shown,  with 
low-amplitude  portions  magnified  by  factors  of  10  and  100. 
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Figure  4  (b).  Narrow-band  pulse,  for  the  case  of  Qo  equal  to  20,  after  it  has  propagated  a 
distance  of  52  times  its  dominant  frequency.  The  pulse  was  constructed  by  convolving  the 
broadband  pulse  with  a  Ricker  wavelet.  The  latter  has  a  dominant  frequency^  equal  to 
twice  the  comer  frequency  of  the  initial  broadband  pulse.  The  initial  portion  is  shown  on 
an  expanded  time  scale. 
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Figure  5.  Q(co)  for  the  numerical  solution  in  the  case  Qo  equal  to  100.  Four  estimates  are 
shown,  corresponding  to  spectral  ratios  at  successive  200  grid-cell  intervals.  The  dotted 

A 

box  delimits  the  usable  computational  bandwidth  and  the  6%  tolerance  interval  about  a  Q 
of  100. 
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Figure  6.  Q{co)  for  the  numerical  solution  in  the  case  Qo  equal  to  20.  Four  estimates  are 
shown,  corresponding  to  spectral  ratios  at  successive  200  grid-cell  intervals.  The  dotted 

box  delimits  the  usable  computational  bandwidth  and  the  6%  tolerance  interval  about  a  Q 
of  19  (which  corresponds  to  a  (2  of  20). 


3  MODEL  FOR  NONLINEAR  WAVE  PROPAGATION  DERIVED  FROM 
ROCK  HYSTERESIS  MEASUREMENTS 

Heming  Xu,  Steven  M.  Day,  and  Jean-Bemard  H.  Minster 
3.1  ABSTRACT 

We  develop  a  method  for  modeling  nonlinear  wave  propagation  in  rock  at 
intermediate  strain  levels,  i.e.,  strain  levels  great  enough  that  nonlinearity  cannot  be 
neglected,  but  low  enough  that  the  rock  does  not  incur  macroscopic  damage.  The 
constitutive  model  is  formulated  using  a  singular-kernel  endochronic  formalism,  and  this 
formulation  is  shown  to  satisfy  a  number  of  general  observational  constraints,  including 
producing  a  power  law  dependence  of  attenuation  (Q'^)  on  strain  amplitude.  Once  the 
elastic  modulus  is  determined,  and  a  second  parameter  fixed  to  give  linear  in  the  strain 
amplitude,  the  model  has  2  remaining  free  parameters.  One  of  these  represents  cubic 
anharmonicity,  and  we  set  in  to  agree  with  laboratory  observations  of  harmonic  distortion. 
The  other  parameter  controls  the  amount  of  hysteresis,  and  it  is  set  to  approximate  stress- 
strain  curves  measured  in  laboratory  uniaxial  stress  experiments  on  Berea  sandstone  by 
Boitnott  and  Haupt.  The  constitutive  equations,  though  fundamentally  nonlinear  and  rate- 
independent,  have  a  superficial,  formal  resemblance  to  viscoelasticity,  which  we  exploit  to 
produce  an  efficient,  stable  numerical  algorithm.  We  solve  ID  wave  propagation  problems 
for  this  constitutive  model  using  both  finite  difference  and  pseudospectral  methods.  These 
methods  are  shown  to  reproduce,  to  high  precision,  analytical  results  for  quasi-harmonic 
wave  propagation  in  a  nonlinearly  elastic  medium. 

Application  of  the  Berea  sandstone  model  to  quasi-harmonic  wave  propagation 
shows  several  departures  from  results  obtained  with  nonlinear  elasticity.  The  Berea  model 
shows  more  rapid  decay  with  distance  of  the  fundamental  frequency  component,  due  to 
nonlinear,  amplitude-dependent  attenuation,  than  does  nonlinear  elasticity.  The  Berea 
model  also  shows  enhanced  excitation  of  the  order  3  harmonic,  in  agreement  with 
laboratory  observations.  In  addition,  the  growth  with  propagation  distance  of  the 
harmonics  of  the  source  excitation  shows  a  saturation,  relative  to  the  nonlinear  elasticity 
results.  This  behavior  reflects  the  competing  effects  of  amplitude  growth  via  energy 
transfer  from  the  source  frequency,  and  energy  dissipation  due  to  hysteresis,  the 
dissipation  increasing  as  the  harmonic  amplitude  grows.  In  additional  numerical 
experiments,  we  find  that  a  two-frequency  source  function  generates  harmonics  with 
frequencies  which  can  be  expressed  as  linear  combinations  of  integer  multiples  of  the  two 
source  frequencies,  in  agreement  with  published  laboratory  results  for  other  solids. 
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Boitnott  and  R.  Haupt  of  New  England  Research,  Inc.  (Boitnott,  1992;  1993).  The  model 
automatically  reduces  to  linear  elasticity  at  low  strain.  The  constitutive  functional  is  recast 
into  a  form  well  adapted  to  large  computations,  in  that  the  model  state  is  specified  by  a 
small  set  of  memory  variables  which  can  be  time-stepped  by  a  simple,  stable  finite 
difference  method.  The  momentum  equation  is  also  time  stepped  by  finite  differencing, 
with  the  spatial  derivatives  approximated  by  a  Fourier  pseudospectral  method. 

We  begin  with  overviews  of  the  general  assumptions  upon  which  the  endochronic 
formalism  rests  and  some  generic  observational  constraints  on  constitutive  models  for  rock. 
Then  we  outline  the  specifics  of  the  endochronic  model,  demonstrate  that  the  formalism  is 
consistent  with  the  generic  observational  constraints,  and  show  the  model  fit  to  Berea 
sandstone  data.  We  then  describe  its  numerical  implementation  for  wave  propagation 
applications  and  show  that  this  method  is  capable  of  giving  accurate  solutions  to  transient 
wave  propagation  problems.  Finally  we  show  that  the  model  predicts  harmonic  distortions 
during  wave  propagation  which  are  comparable  to  those  observed  in  laboratory  wave 
propagation  experiments. 

3.3  BACKGROUND 

Constitutive  Functionals 

We  fit  laboratory  data  and  simulate  wave  propagation  using  an  elastic-plastic 
constitutive  model  of  the  so-called  endochronic  type.  We  begin  this  section  with  a  brief 
overview  of  constitutive  models,  in  order  to  indicate  the  general  assumptions  which 
underlie  the  endochronic  formalism. 

We  restrict  consideration  to  so-called  simple  materials,  i.e.,  materials  for  which  the 
stress-strain  relationship  is  local,  in  the  sense  that  stress  at  material  point  X  is  a  functional 
of  the  deformation  gradient  only  (and  not,  for  example,  dependent  on  higher  spatial 
derivatives  of  the  displacements).  Assuming  the  principle  of  frame  indifference  applies,  the 
second  Piola-Kirchhof  stress  tensor  at  time  t,  n(0,  is  then  a  functional  of  the  (Lagrangian) 
strain  tensor  time  history,  E(r): 

7c(X,0  =  F[E(X,O,— (1) 

When  the  displacement  gradients  are  small,  we  can  use  the  Cauchy  stress  tensor  S 
in  place  of  K  and  the  infinitesimal  strain  tensor  6  in  place  of  E, 

S(X.0  =  F[£(X,r'),— <?'<?].  (2) 

This  approximate  formulation  still  permits  nonlinearity  of  the  constitutive  functional,  but 
neglects  geometrical  nonlinearities.  For  most  applications  to  earthquake  and  explosion 
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ground  motion,  this  approximation  is  adequate,  and  we  adhere  to  it  in  the  current  study. 
Important  exceptions  would  be  motions  very  near  the  source,  as  well  as  ground-motion- 
induced  soil  liquefaction,  both  of  which  are  beyond  our  scope. 

If  memory  effects  are  suppressed  in  (1)  or  (2),  so  that  stress  is  simply  a  function  of 
the  strain,  we  have  the  special  case  of  nonlinear  elasticity.  In  this  case,  there  is  no 
hysteresis,  and  wave  propagation  entails  no  loss  of  energy,  but  there  will  in  general  be 
transfer  of  energy  among  frequency  components  of  the  wavefield. 

If  memory  effects  are  retained,  but  the  functional  in  (2)  is  time  shift  invariant  and 
linear,  we  have  the  special  case  of  linear  viscoelasticity,  (2)  reducing  to  a  convolution 

S(r)  =  M(0*^.  (3) 

at 


where  M  is  the  (fourth-order  tensor  valued)  time-dependent  modulus.  The  material  will 
exhibit  hysteresis,  with  elliptical  loops  under  monochromatic  loading.  Superposition 
applies,  and  there  is  no  energy  transfer  among  frequency  components.  If  memory  effects 
are  also  suppressed  in  (3),  we  have  linear  elasticity. 

We  can  rewrite  the  general  functional  (1)  in  the  form 

it(X,0  =  f[E(X,^');i(X,|');-~  S  «'  <  |(X,r)]  (4) 

where  ^  is  the  strain  path  length,  the  increment  of  which  is  given,  in  terms  of  positive 
definite  tensor  g,  by 


4=dE:g:dE.  (5) 

If  the  stress  is  a  functional  of  strain,  expressed  as  a  function  of  but  has  no  explicit 
dependence  on  the  strain  rate  the  rheology  is  said  to  be  rate-independent,  and  (4) 
reduces  to 


jc(X,r)  =  F[E(X,|0;—  ^  ^  l(X,r)].  (6) 

In  this  case,  the  wavefield  (in  a  homogeneous  medium)  obeys  a  simple,  well-known 
scaling  relationship:  if  ui(X,t)  is  a  displacement  field  obeying  the  equations  of  motion, 
then  a  positive  parameter  X  generates  a  family  of  solutions  ux(X,0  given  by 

u^(X,t)  =  Xu,iX-^XX't).  (7) 
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The  scaling  relationship  remains  valid  whether  or  not  geometrical  nonlinearities  are 
significant.  Linear  and  nonlinear  elasticity  and  classical  plasticity  are  examples  of  rate- 
independent  rheologies,  but  linear  viscoelasticity  is  not  (apart  from  the  special  case  of  linear 
elasticity).  Rate-independent  materials  may  exhibit  hysteresis,  and  the  wavefield  will  lose 
energy  during  wave  propagation.  In  contrast  to  the  viscoelastic  case,  however,  the 
hysteresis  loops  are,  in  general,  not  elliptical. 

Here  we  adopt  the  so-called  endochronic  plasticity  formalism  of  Valanis  and  Read 
(1982),  which  employs  a  special  case  of  the  rate-independent  functional  (6).  This 
formalism  partitions  strain  increments  into  elastic  and  inelastic  ("plastic")  components,  such 
that 


rf7t  =  G(tflE-dEa  (8) 

where  G  is  an  elastic  tensor  and  dEP  denotes  the  plastic  strain  increment.  A  constitutive 
functional  is  constructed  in  which  the  strain  in  (6)  is  parameterized  in  terms  of  some  plastic 
strain  path  length  z 


7i(X,0  =  F[E(X,z(^');0  <^'<  ^(X,0].  (9) 

The  functional  (9)  is  then  written  in  the  form  of  an  integral  of  a  shift-invariant  (in  z)  kernel 
over  the  plastic  strain  path.  The  theory  thus  has  some  formal  resemblance  to  viscoelasticity, 
but  with  the  role  of  time  as  the  integration  variable  in  the  constitutive  functional  replaced  by 
z  (as  a  result  of  this  formal  similarity,  z  is  sometimes  referred  to  as  an  "intrinsic  time," 
which  is  responsible  for  the  designation  "endochronic").  The  formalism  is  inherently  rate- 
independent,  the  constitutive  functional  having  no  explicit  dependence  upon  the 
deformation  rate. 

We  next  outline  some  general  observational  constraints  that  apply  to  constitutive 
models  for  rock  at  intermediate  strain,  then  present  a  suitable  one-dimensional,  small-strain 
implementation  of  (8-9).  The  resulting  constitutive  model  can  be  shown  to  have  a  number 
of  features  consistent  with  the  observational  constraints  outlined  below,  though  it  is  by  no 
means  unique  in  that  respect.  We  can  also  exploit  the  superficial  resemblance  to 
viscoelasticity  to  obtain  a  concise,  efficient,  and  stable  numerical  formulation,  rendering  the 
model  very  suitable  for  wave  propagation  simulations. 

Observations  of  Nonlinear  Attenuation 

At  low  strain  levels,  less  than  about  10'^,  viscoelasticity  has  been  found  to  provide 
an  adequate  framework  for  interpreting  energy  losses  in  seismic  wave  propagation, 
including  the  frequency  dependence  of  those  losses  (e.g..  Minster,  1980).  At  intermediate 
strain  levels,  however,  the  situation  is  different.  In  this  strain  regime,  the  inadequacy  of 
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linear  viscoelasticity  can  be  seen  in  both  amplitude  dependence  and  frequency  scaling  of  the 
energy  loss  mechanism.  While  the  deformation  is  not  sufficient  to  cause  macroscopic 
failure,  it  is  large  enough  that  preexisting  microcracks  may  propagate,  new  ones  may  be 
created  under  stress,  and  even  microplasticity  may  be  induced  (Mashinsky,  1994). 
Because  the  rheology  depends  on  amplitude,  wave  propagation  is  characterized  by 
nonlinearity  and  reflects  the  influence  of  inelastic  mechanisms.  In  the  intermediate  strain 
regime,  these  mechanisms  dominate  and  control  the  rate  and  mode  of  energy  dissipation, 
but  in  the  low  strain  regime,  their  effects  are  masked  by  linear,  amplitude-independent 
viscoelastic  mechanisms.  Since  high  to  intermediate-level  strains  occur  near  seismic 
sources,  the  proper  representation  of  nonlinear  constitutive  equations  for  rocks  in  this 
regime  is  a  potentially  important  ingredient  of  quantitative  source  models. 

Seismic  wave  attenuation  at  intermediate  strain  levels  has  been  observed  to  be 
amplitude  dependent  on  the  basis  of  both  laboratory  and  field  experiments.  Apparent 
as  estimated  from  resonance  bar  and  pulse  transmission  experiments  in  the  laboratory,  has 
been  found  to  be  approximately  proportional  to  strain  amplitude  for  strains  exceeding  about 
10^  (e.g.,  Winkler  et  al.,  1979;  Mavko,  1979;  Johnston  and  Toksoz,  1980;  Stewart  et  al., 
1983;  Bonner  and  Wanamaker,  1990).  From  analysis  of  seismic  data  from  field 
experiments  with  explosive  sources.  Minster  and  Day  (1986)  concluded  that  propagation 
losses  in  salt  are  amplitude  dependent  at  intermediate  strain  levels,  and  that  the  apparent  Q~^ 
is  compatible  with  the  linear  dependence  of  apparent  Q~^  on  strain  amplitude  commonly 
observed  in  laboratory  experiments. 

When  propagation  losses  in  the  intermediate  strain  range  are  analyzed  in  the  context 
of  linear  viscoelasticity,  nonlinearity  can  induce  an  apparent  frequency  dependence  of  Q. 
For  example,  Wortman  and  McCartor  (1987)  estimated  Q(f)  from  spectral  ratios  of  ground 
motion  recordings  (maximum  strain  ~4  x  10‘4)  from  the  Salmon  explosion  (5.3  kt)  in  salt. 
They  found  Qif)  decreasing  sharply  with  decreasing  frequency  below  the  6  Hz  comer 
frequency  of  the  explosion.  They  performed  the  same  analysis  for  one  of  the  much  smaller 
Cowboy  explosions  (corner  frequency  ~100  Hz),  also  in  salt,  with  recordings  at  the  same 
scaled  range  and  similar  maximum  strain  level.  In  the  latter  case,  Wortman  and  McCartor 
found  an  absorption  band  of  very  similar  shape,  but  shifted  up  in  frequency  in  proportion 
to  the  corner  frequency  ratio  of  the  sources.  This  result  is  incompatible  with  the 
assumption  of  linear  viscoelasticity.  However,  the  shift  in  the  apparent  absorption  band  is 
precisely  the  scaling  that  would  be  predicted  for  a  nonlinear,  rate-independent  rheology 
(i.e.,  Eqn  7). 

Observations  of  Harmonic  Distortion 

The  amplitude  dependence  of  attenuation  at  intermediate  strains  is  most  likely 
associated  with  friction  along  microcracks  and  joints  (Boitnott,  1992;  1993).  Considerable 
experimental  and  theoretical  interest  is  developing  in  this  and  other  nonlinear  effects 
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associated  with  the  propagation  of  sound  waves  in  nonlinear  solids.  Continuous  acoustic 
wave  phase  detections  of  wave  propagation  in  solids,  including  rocks,  have  found  that 
when  a  sinusoidal  ultrasonic  wave  of  a  given  frequency  and  of  sufficient  amplitude  is 
introduced,  the  fundamental  wave  distorts  as  it  propagates  and  the  higher  harmonics  of  the 
fundamental  frequency  are  generated  (Breazeale  and  Ford,  1965;  Hirata  et  al.,  1965; 
Thompson  et  al.,  1976;  Thompson  and  Tiersten,  1977;  Meegan  et  al.,  1993).  In  summary, 
ultrasonic  experiments  have  identified  these  characteristics  of  wave  propagation  and 
interaction:  1)  harmonics  of  all  orders  are  generated  for  a  single  frequency  sinusoidal 
source;  2)  the  fundamental  component  decays  away  from  source  and  the  second  harmonic 
amplitude  increases  linearly  with  the  propagation  distance  and  quadratically  with  both 
source  amplitude  and  the  source  frequency  (e.g.,  Meegan  et  al.,  1993  for  Berea  sandstone 
samples).  3)  two-frequency  sound  beams  are  found  to  interact  in  nonlinear  solids  so  as  to 
produce  frequencies  of  combinations  (sums  and  differences  of  integer  multiples)  of  the  two 
driving  frequencies  (Bonner  and  Wanamaker,  1990,  1991;  Johnson  et  al,  1987,  1991; 
Johnson  and  Shankland,  1989).  Since  energy  transfer  between  frequency  bands  can  not 
occur  in  a  linear  system,  such  observations  are  clearly  diagnostic  of  nonlinear  rheology  at 
intermediate  strain  levels.  The  potential  use  of  these  nonlinear  effects  as  a  means  of 
exploring  the  early  stages  of  the  plastic  deformation  and  mechanism  of  dislocation  damping 
has  been  demonstrated  in  the  experiments  of  Hirata  et  al.  (1965)  and  Morris  et  al.  (1979). 

The  harmonic  distortions  during  wave  propagation  may  result  from  at  least  two 
causes.  The  lattice  itself  can  be  anharmonic,  or  dislocations  vidthin  the  solid  can  produce 
such  nonlinear  effects  (Breazeale  and  Ford,  1965).  The  formalism  of  nonlinear  elasticity 
has  been  applied  to  describe  the  consequences  of  lattice  anharmonicity  (Breazeale  and  Ford, 
1965),  which  implies  material  imperfection  on  the  atomic  scale.  Nonlinear  effects  due  to 
dislocation  displacements  have  been  considered  as  well  (Hirata  et  al.,  1965).  Nonlinear 
interactions  of  frequency  components  in  large-amplitude  elastic  waves  in  rocks  were 
analyzed  with  perturbation  expansions  and  Green  function  techniques  in  terms  of  nonlinear 
elasticity  (McCall,  1994;  Thompson  and  Tiersten,  1977).  For  a  monochromatic 
longitudinal  plane  wave,  for  example,  the  excitation  of  higher  harmonics  can  be  shown  to 
scale  with  a  dimensionless  coefficient  fi,  which  represents  the  cubic  anharmonicity,  and  is 
defined  in  terms  of  the  second  and  third  order  elastic  moduli.  For  example,  when  the 
amplitudes  of  higher  harmonics  measured  in  wave  propagation  experiments  in  Berea 
sandstone  were  interpreted  in  the  context  of  nonlinear  elasticity  theory,  the  coefficient  p 
was  estimated  to  be  about  7x103  (Meegan  et  al.  1993).  Other  models  have  incorporated 
asymmetries  between  loading  and  unloading  moduli  to  represent  solids  with 
microplasticity;  for  example,  Nazarov  et  al.  (1988)  propose  such  a  model  to  account  for  the 
experimental  results  on  harmonics  generation  by  standing  waves  on  a  copper  rod. 
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Observations  of  Stress-Strain  Curves 

Laboratory  stress-strain  curves  under  cyclic  loading  characteristically  demonstrate 
the  following  features  which  a  successful  model  must  emulate:  1)  Hysteresis  occurs, 
implying  energy  loss,  and  the  effective  Q  characterizing  the  dissipation  is  strain-amplitude 
dependent.  2)  Hysteresis  loops  are  cusped  at  reversal  points,  rather  that  elliptical  (  as 
would  typify  linear  anelastic  behavior);  3)  No  yield  surface  is  evident  in  the  loading  curves, 
at  least  for  strains  up  to  10  '^.  4)  Upon  reversal  of  strain  path,  the  tangent  modulus  is 
roughly  equal  to  the  instantaneous  elastic  modulus  (  Gordon  and  Davis  ,1968;  Liu  and 
Peselnick,  1979;  McKavanagh  and  Stacey,  1974;  Minster  et  al.,  1991).  These,  hysteresis 
loops  indicate  that  weak  microplasticity  is  created  during  the  cyclic  loading  processes  at  the 
intermediate  strain  levels.  Attempts  to  simulate  this  kind  of  hysteretic  behavior  using  linear 
anelasticity  as  the  point  of  departure  fail  to  replicate  one  or  more  of  the  foregoing  features 
(Minster  et  al.,  1991;  Day  et  al.,  1992).  The  observed  behavior  can  be  represented  in  terms 
of  an  ensemble  of  many  simple  hysteretic  units  (McCall  and  Guyer,  1994;  Guyer  et  al., 
1995;  1997),  and  the  properties  of  the  ensemble  may  provide  insight  into  the 
micromechanics  underlying  the  hysteresis,  but  such  representations  do  not  of  themselves 
provide  a  formulation  suitable  for  large  numerical  simulations  of  wave  propagation. 

3.4  IMPLEMENTATION  OF  THE  ENDOCHRONIC  MODEL 
Singular-Kernel  Formulation 

We  assume,  following  Valanis  and  Read  (1982),  that  (9)  can  provide  an  adequate 
phenomenological  description  of  intermediate-strain  behavior  of  rock  when  restricted  to 
shift-invariant  integral  form.  Here  we  will  consider  only  scalar  waves  in  one  dimension,  in 
the  small  displacement  gradient  approximation  (Heming  et  al.,  1998,  in  preparation, 
describe  the  generalization  to  3D).  Denoting  the  single  stress  and  strain  components  by  S 
and  e,  respectively,  (8)  and  (9)  become 

dS  =  Gide-de'’)  ,  (10) 


and 


(11) 


where  G  is  the  elastic  modulus,  and  the  intrinsic  time  increment  dz  is  defined  in  terms  of 
the  plastic  strain  increment  rffP  by 
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(12) 


If  the  kernel  function  is  given,  then,  given  an  increment  of  strain,  the  stress  increment  can 
be  calculated  from  the  Equations  10-12,  as  we  will  show  below. 

If  the  kernel  function  is  weakly  singular  at  z=0,  taking  the  form 

A:(z)~z-“  (13) 

where  a  is  positive,  then  the  constitutive  functional  satisfies  the  4  observational  constraints 
on  hysteresis  loops  noted  in  the  previous  section.  In  particular,  Valanis  and  Read  (1982) 
show  that  (13)  implies  that  both  initial  loading  and  subsequent  unloading  and  reloading  at 
stress  reversals  will  occur  along  the  elastic  slope  G  in  the  stress-strain  plane,  in  general 
agreement  with  the  behavior  of  rocks  at  intermediate  strain  levels. 

Amplitude  Scaling  of  Q, 

We  can  also  demonstrate  that  the  singular  kernel  ensures  power  law  dependence  of 
Q-^  on  strain  amplitude  at  intermediate  strain  levels,  in  agreement  with  the  experimental 
observations  reviewed  earlier.  We  consider  a  cycle  through  a  particular  hysteresis  loop 
Si(z).  Construct  Si(z)  by  the  following  sequence  of  steps.  First  load  monotonically  until 
the  stress  reaches  some  prescribed  maximum  value  Sim-  At  that  point  z  equals  some  value 
zi,  so  Sim=Si(zi).  Then  unload  and  load  in  the  reverse  direction  to  -Sim,  corresponding  to 
some  z  =Z2,  and  finally  reverse  again  and  reload  back  to  Sim>  corresponding  now  to  some 
z=Z3.  The  inelastic  work  Wi  done  (per  unit  volume)  by  the  load  over  the  loop  Si  is 

W,=|s,(z)^dz  .  (14) 

Consider  a  similar  cycle,  but  with  the  reversal-point  stress  amplitudes  scaled  such  that  z  at 
each  reversal  point  changes  by  some  positive  scale  factor  y,  to  give  reversals  at  yzi,  i=l,  2, 
3.  It  is  clear  from  (11-13)  that  the  new  stress  function,  denoted  is  a  scaled  version 
of  Si, 


S,(z)  =  y‘-“S;(y-‘z), 

By  (14)  and  (15),  the  inelastic  work  Wyover  the  scaled  loop  is 


(15) 


(16) 
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We  define  Q'^  for  a  loop  of  a  given  amplitude,  by  analogy  with  viscoelastic  Q'^,  as  l/2jt 
times  the  ratio  of  net  inelastic  work  done  over  the  cycle  to  average  strain  energy  during  the 
cycle.  We  assume  that  the  inelastic  work  is  small  compared  with  the  average  strain  energy, 
i.e.,  Q.»\,  in  which  case  the  peak  strain  is  approximately  proportional  to  peak  stress,  and 
the  average  strain  energy  is  approximately  proportional  to  the  product  of  peak  strain  and 
peak  stress.  Then,  combining  (15)  and  (16),  we  find 

(2,-'  =  yV.  (17) 

where  Q\  denotes  Q  for  the  reference  loop  5i.  Using  (15)  we  can  rewrite  the  scale  factor 
in  (17)  in  terms  of  maximum  stress  amplitude,  and  with  the  low  loss  approximation  noted 
above,  we  can  replace  maximum  stress  amplitude  with  maximum  strain  amplitude. 
Denoting  the  maximum  strain  amplitude  by  em.  and  value  of  £m  for  the  reference  loop  by 
eim,  we  obtain 


/p  \°A-. 


(18) 


where  we  have  now  written  0(eim)  to  denote  Q\.  Thus,  (2*1  has  approximately  power  law 
dependence  on  strain  amplitude,  with  scaling  exponent  a/(l-a).  For  the  choice  a  =  1/2, 
we  have  an  approximately  linear  dependence  of  on  strain  amplitude  em>  in  agreement 
with  numerous  laboratory  observations. 

Memory  Variable  Approximation 

The  integral  form  for  the  stress,  (1 1),  is  intractable  for  numerical  wave  propagation 
computations,  since  the  integral  is  performed  over  the  entire  strain  history.  This  limitation 
is  analogous  to  that  faced  in  the  numerical  implementation  of  viscoelasticity  (Day  and 
Minster,  1984),  and  our  solution  is  similar.  That  is,  recognizing  the  convolutional  form  of 
(11),  we  can  replace  the  integral  with  an  approximation  to  the  corresponding  differential 
operator.  This  operator  follows  once  the  kernel  function  is  expanded  in  terms  of  decaying 
exponentials  in  z  (e.g.,  Valanis  and  Read  1982).  We  develop  such  an  expansion  for  the 
kernel  function  C/z^  using  the  method  of  Prony  (e.g.,  Marple,  1987),  leading  to  an 
approximation  of  the  form 


/(:(z)  =  Xv“'‘ 

r=l 


(19) 
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Once  K{z)  is  approximated  by  the  exponential  sum  (19),  (11)  gives  stress  in  the 

form 


where 


dz' 


(20) 


(21) 


The  p  constants  ax  and  Ax,  are,  respectively,  the  exponents  and  coefficients  from  the  order  p 
Prony  expansion  of  K.  Then  (21)  may  be  differentiated  to  yield  the  following  linear,  first- 
order  differential  equations  for  the  p  internal,  or  "memory,"  variables  qx'. 


dz 


'-  +  a,q. 


r=l,2,...p. 


(22) 


The  integral  form  (11)  has  been  changed  to  incremental  form  (22),  and  the  memory 
requirement  reduced  to  storing  the  p  memory  variables.  It  is  no  longer  necessary  to  store 
the  full  strain  history,  nor  even  to  explicitly  save  the  reversal  points.  In  practice,  we  have 
found  that  3  to  4  memory  variables  per  stress  node  provides  an  approximation  adequate  for 
simulation  of  wave  propagation  in  the  intermediate-strain  range. 

Numerical  Solution 


The  constitutive  equations  now  consist  of  (22),  (10)  and  (12),  which  may  be  solved 
numerically.  A  convenient  numerical  scheme  for  the  solution  of  such  a  system  is  given  by 
Murakami  and  Read  (1989).  We  found  it  necessary  to  generalize  the  Murakami  and  Read 
approach  somewhat,  retaining  higher  order  accuracy  in  the  time  stepping  of  (22).  Taking 
de^Jdz  to  be  constant  &'’/&  ,  integrating  (21)  from  zq  to  zo+5z ,  we  obtain 


7,(Zo  +  Sz)  =  q,(Zo)e  +  4^(1  - 


Combining  (23)  with  (10),  (12)  and  (20)  yields 


GSz  + 


A,  Se" 

r=lV«r 


Se” 

Sz 


(23) 


(24) 
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Isolating  5fP  and  squaring  leads  to  an  equation  for  5z  which  is  easily  solved  numerically. 
Then  (23)  can  be  stepped  forward  by  5z  to  update  the  memory  variables,  after  which  the 
updated  stress  is  evaluated  from  (20). 

3.5  FITTING  QUASI-STATIC  LABORATORY  DATA 

We  have  shown  that  the  singular-kernel  model  is  consistent  with  the  main  generic 
features  of  rock  nonlinearity.  Nonetheless,  the  model  is  essentially  phenomenological,  in 
that  it  is  not  derived  from  micromechanical  considerations.  The  constant  factor  C  which 
scales  the  kernel  function  must  be  calibrated  by  matching  the  model  to  laboratory  stress- 
strain  observations.  Furthermore,  the  range  of  z  over  which  the  Prony  expansion  of  the 
kernel  is  optimized  will  affect  the  shape  of  hysteresis  loops  at  larger  strain  amplitudes, 
which  introduces  a  further  degree  of  freedom  which  may  be  constrained  by  fitting 
laboratory  data.  Here  we  fit  the  model  to  stress-strain  data  from  quasi-static  uniaxial  tests 
of  Berea  sandstone  (G.  Boitnott  and  R.  Haupt,  New  England  Research,  personal 
communication;  Boitnott,  1992;  1993). 

Typical  raw  laboratory  data  in  the  form  of  stress  and  strain  histories  are  often  rather 
noisy  at  the  relatively  low  strain  levels  of  interest  to  us,  and  require  filtering.  Simple  low- 
pass  filtering  smoothes  the  cusps,  thereby  masking  the  onset  of  nonlinear  behavior,  and 
affecting  the  measurement  of  moduli  at  and  near  the  cusps.  We  have  therefore  developed  a 
technique  to  filter  separately  the  loading  and  unloading  portions  of  the  loops.  It  relies  on 
the  construction  of  a  longer  time  series  out  of  a  half-loop — that  is  a  portion  of  stress-strain 
history  between  two  reversals,  in  which  both  stress  and  strain  are  monotonic.  This  is  done 
by  extending  it  in  both  directions  with  versions  of  itself,  rotated  by  ±7C  about  its  end 
points.  The  extended  time  series  is  then  de-meaned,  de-trended,  and  low-pass  filtered 
using  a  zero-phase  filter  to  avoid  introduction  of  a  phase  shift.  The  filtered  version  is 
truncated  to  the  original  length  after  restoring  trend  and  mean,  and  the  stress-strain  path 
reconstructed  by  concatenation  of  filtered  segments.  The  rotation  by  ±7t  of  the  extensions 
has  the  advantage  of  preserving  the  continuity  of  the  time  series  and  its  derivative,  thereby 
limiting  undesirable  end  effects.  It  is  important  to  avoid  introducing  cusps  into  hysteresis 
loops  when  they  are  not  present,  and  to  avoid  smoothing  through  a  cusp  or  changing  its 
angle,  when  one  is  present.  Too  strong  a  filter  will  change  the  slope  near  the  end  points, 
and  therefore  introduce  a  fictitious  cusp.  Misapplication  of  the  technique  is  detectable 
because  this  will  create  overlaps  or  gaps  between  successive  segments.  The  technique  gives 
very  satisfactory  results  for  noise  levels  as  large  as  10  percent,  as  we  have  been  able  to 
verify  using  synthetic  loops  contaminated  with  additive  noise.  This  approach  facilitates 
considerably  the  estimation  of  the  tangent  modulus,  particularly  near  the  loop  ends  where 
noise  contamination  is  most  worrisome. 
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We  first  find  the  two  constants  G  and  C.  As  shown  by  Valanis  and  Read  (1982), 
if  the  kernel  function  takes  the  form  Ar(z)=  C  z‘“*  then  during  the  initial  loading  process  the 
stress  and  strain  have  an  explicit  relationship  which  is 


(25) 


We  use  a  least  squares  fit  to  stress  and  strain  values  at  hysteresis  loop  reversal  points,  for 
varying  maximum  stress  levels,  to  determine  the  constants  G  and  C.  This  is  shown,  for 
example,  in  Fig.  la  for  the  Berea  sandstone  data  of  Boitnott  and  Haupt.  We  then  use  the 
Prony  expansion  to  get  and  .  A  comparison  between  the  exact  and  approximate  kernel 
functions  for  a=l/2,  using  4  terms  in  the  expansion,  is  shown  in  Fig.  lb.  The  exponent 
is  fixed  at  a=l/2  since  that  value  yields  a  linear  dependence  of  Q'^  on  strain  amplitude,  as 
shown  by  (18),  a  result  consistent  with  experimental  data.  Table  1  lists  all  model 
parameters  resulting  from  the  fitting.  Fig.  Ic  and  Id  show  a  series  of  calculated  hysteresis 
loops  for  the  exact  and  approximate  kernels,  respectively.  There  is  a  saturation  effect 
associated  with  the  Prony  approximation  when  the  strain  level  is  sufficiently  high  (>~10'^) 
which  is  not  evident  in  these  figures.  This  is  due  primarily  to  the  exponential  decay  of  the 
Prony  approximation  at  large  values  of  zjc,  compared  with  the  power  law  decay  of  the  exact 
kernel.  Further  terms  in  the  Prony  approximation  (and  therefore  additional  memory 
variables)  would  be  required  to  maintain  the  power  law  decay  of  the  exact  kernel  to  higher 
strain  levels.  Thus,  data  collected  at  very  large  strains  may  require  a  higher  order 
approximation  to  achieve  a  good  fit. 

We  test  the  resulting  model  by  computing  synthetic  hysteresis  loops  for  a  range  of 
maximum  strain  levels,  which  we  compare  with  the  corresponding  laboratory  curves.  The 
comparison  for  3  different  peak  strain  levels  is  shown  in  Fig.  2.  As  illustrated  in  Figure  2, 
this  model  is  very  successful  for  modeling  high  signal-to-noise  data  sets  such  as  the  Berea 
sandstone  loops  for  strain  levels  between  7-40x10'^.  The  theoretical  loops  simulate  the 
amplitude  dependence  of  Q  and  the  non-elliptical  loop  shapes,  including  cusps  at  the  ends. 

There  is  no  reason  that  the  elastic  modulus  G  in  the  endochronic  model  need  be 
constant.  The  modulus  could  be  made  strain  dependent,  thereby  adding  a  component  of 
nonlinear  elasticity  to  the  Berea  model.  Doing  so  does  not  discemibly  improve  the  fit  to  the 
quasi-static  data,  however.  Our  preliminary  model  therefore  uses  a  strain-independent 
modulus.  But  when  we  apply  the  preliminary  model  to  wave  propagation  in  the  next 
section,  comparison  of  the  simulations  with  laboratory  narrowband  wave  propagation  data 
reveals  a  deficiency  in  the  simulation  of  even  harmonics  of  the  source  frequency.  This 
deficiency  is  readily  corrected  by  giving  the  elastic  modulus  a  weak  strain  dependence.  The 
generation  of  harmonics  thus  provides  a  particularly  sensitive  further  test  of  the  model. 

In  its  ability  to  match  observations  of  rock  hysteresis,  this  approach  represents  a 
substantial  improvement  over  attempts  to  simulate  amplitude-dependent  attenuation  using 


45 


variants  of  viscoelasticity  (Minster  and  Day,  1986;  Minster  et  al.,  1991;  Day  et  al.,  1992). 
Although  purely  phenomenological,  the  endochronic  approach  has  the  decided  advantage 
that  it  readily  reduces  to  a  relatively  simple  differential  equations  which  are  easily 
programmed  and  solved  numerically.  The  algorithm  described  above  is  sufficiently  stable 
and  efficient  that  it  can  be  applied  to  compute  stress-strain  behavior  in  large  numerical  wave 
propagation  simulations.  Efficiency  and  stability  of  the  algorithm  are  critical  in  large  2D  or 
3D  simulations,  in  which  10^^  or  more  evaluations  of  the  constitutive  functional  will 
typically  be  necessary.  We  demonstrate  the  method  for  the  ID  case. 

3.6  NONLINEAR  WAVE  PROPAGATION 
Methods 

Standard  Fourier  synthesis  methods  for  simulating  ground  motion  fail  when  rock  or 
soil  nonlinearity  is  introduced  into  the  continuum  model.  While  considerable  insight  into 
the  resultant  wave  propagation  phenomena  has  been  obtained  through  approximations 
(e.g.,  McCall,  1994),  general  treatment  of  the  problem  requires  volume  discretization 
methods  such  as  finite  differencing.  In  the  intermediate  strain  regime,  nonlinear  losses, 
pulse  distortions  and  harmonic  generation  have  been  well  documented  for  solids  in  the 
laboratory.  In  this  study,  two  methods  have  been  developed  for  one-dimensional  nonlinear 
wave  propagation  with  the  endochronic  model:  (i)  a  one  dimensional  fourth-order  finite 
difference  method  which  is  fourth-order  accurate  in  space  (second  order  in  time),  and  (ii)  a 
one-dimensional  pseudospectral  method,  with  second-order  finite  differencing  in  time 
(e.g.,  Kosloff  and  Baysal,  1982;  Fornberg,  1987;  Daudt  et  al.,  1989).  In  the  latter 
method,  spatial  derivatives  are  approximated  to  infinite  order  via  discrete  Fourier 
transforms.  We  have  found  that  these  methods,  both  of  which  use  spatial  differencing  of 
higher  than  second  order  accuracy,  work  well  in  this  application,  just  as  they  known  to  in 
linearly  elastodynamic  computations.  This  is  in  part  because  we  have  confined 
consideration  to  sources  with  dominant  wavenumber  several  octaves  below  the  Nyquist 
wavenumber  of  the  grid,  and  restricted  the  computations  to  relatively  low  strain  levels.  As  a 
result,  nonlinear  effects  are  not  strong  enough  to  transfer  appreciable  energy  to 
wavenumbers  approaching  the  Nyquist.  Since  we  are  studying  waves  in  the  intermediate 
strain  regime,  shock  and  macroscopic  failure  phenomena  are  not  important  factors,  and 
losses  from  hysteresis  further  smooth  the  wavefield.  In  this  situation,  the  preponderance  of 
the  computational  load  arises  firom  repeated  evaluations  of  the  constitutive  functional.  The 
computational  advantage  of  using  differencing  of  higher  than  second  order  thus  becomes 
very  attractive,  because  the  number  of  stress  evaluations  for  a  given  level  of  accuracy  can 
be  significantly  reduced;  in  the  case  of  fourth-order  differencing,  for  example,  stress 
evaluations  can  be  reduced  by  approximately  a  factor  of  2n,  relative  to  second-order 
differencing,  where  n  is  the  number  of  spatial  dimensions.  The  pseudospectral  method 
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further  reduces  the  relative  contribution  of  the  constitutive  equations  to  the  total 
computational  load. 

A  further  computational  advantage  of  the  fourth-order  and  pseudospectral  methods 
over  the  second-order  methods  which  may  be  more  familiar  in  nonlinear  computational 
studies  is  the  reduction  of  numerical  dispersion.  Dispersion  associated  with  low-order 
differencing  mainly  affects  high  frequencies,  but  in  nonlinear  computations  (and  in  contrast 
to  linear  elastodynamic  computations),  numerical  artifacts  at  high  frequency  can  couple 
back  into  low  frequency  components  of  the  wavefield,  and  vice  versa.  For  this  reason, 
artificial  viscosity  is  essential  in  second-order  nonlinear  computations  to  continuously 
suppress  the  high  frequency  content.  We  have  found,  for  the  relatively  low  strain  levels 
considered  in  this  study,  that  both  the  fourth-order  and  pseudospectral  methods  obviate  the 
need  for  artificial  viscosity  in  intermediate-strain  computations  with  the  endochronic  model. 
This  simplifies  the  task  of  modeling  and  analyzing  the  spectral  distortions  introduced  by  the 
nonlinear  constitutive  model. 

One  of  the  better  understood  consequences  of  material  nonlinearity  on  waveforms 
is  harmonic  distortion  (e.g.,  Hirata  et  al.,  1965;  Breazeale  and  Ford,  1965).  A 
monochromatic  source  generates  harmonics  which  grow  both  with  propagation  distance 
and  with  source  amplitude.  Response  to  quasi-monochromatic  excitation  is  of  interest  in  its 
own  right  ( for  example,  in  the  field  of  non-destructive  testing  (Bonner  and  Wanamaker, 
1990,1991;  Morris  et  al,  1979));  it  is  also  a  starting  point  for  building  a  good 
understanding  of  spectral  distortions  in  nonlinear  solids  in  general.  Furthermore,  previous 
analytical  work  and  experimental  results  provide  a  starting  point  for  validation  and 
interpretation  of  numerical  computations  (Beresnev  and  Nikolaev,  1988;  McCall,  1994; 
Meegan  et  al,  1993;  Thompson  and  Tiersten,  1977).  For  these  reasons,  we  will  apply  the 
numerical  method  to  quasi-monochromatic  sources  in  the  following  sections. 

We  have  tested  the  two  one  dimensional  wave  propagation  codes  for  the  case  of 
linear  elasticity.  With  both  methods,  we  calculated  the  propagation  of  a  narrowband 
wavelet  consisting  of  a  sinusoid  modulated  with  a  Gaussian  envelope  and  verified  that  the 
fundamental  pulse  remains  invariant  with  distance  and  that  no  higher  harmonics  are 
generated  in  the  spectra. 

Nonlinear  Elasticity 

As  a  verification  of  the  accuracy  of  the  wave  propagation  calculations  in  a  nonlinear 
problem,  we  present  numerical  results  from  nonlinear  elastic  wave  propagation 
calculations.  Much  theoretical  work  has  been  focusing  on  wave  propagation  in  nonlinear 
elastic  materials  (e.g.,  McCall,  1994).  We  compare  the  numerical  solutions  from  finite 
difference  and  pseudospectral  methods  with  the  approximate  analytical  solutions  obtained 
using  McCall’s  perturbation  expansion  method.  For  this  comparison,  the  stress-strain 
relationship  is 


47 


a  =  G(e)e , 


(26) 


where 


C(£)  =  Go(l  +  /fe),  (27) 

and  P  is  the  cubic  anharmonicity  coefficient.  We  have  set  P  to  5000,  which  is  in  the  range 
of  10^-10^  for  rock  suggested  by  Johnson  and  McCall  (1994),  and  near  the  estimate  of 
Meegan  et  al.  (1993)  for  Berea  sandstone.  The  source  is  taken  to  be  a  longitudinal,  planar 
force  distribution  (i.e.,  point  force  in  ID).  For  comparison  with  analytical  results  for 
monochromatic  excitation,  we  use  a  narrowband  source  time  function,  a  Gaussian- 
modulated  sinusoid 


Fit)  =  Aexp[-(r  -  2wTf  /iwTf]smi2m/T),  (28) 

where  T  is  the  dominant  period  and  w  is  a  measure  of  the  duration  of  the  source,  in  units  of 
the  dominant  period.  Eqn.  28  has  a  Gaussian-shaped  spectral  peak  with  amplitude  equal  to 
the  time  domain  peak  amplitude  times  the  factor  4itwTf2.  To  ensure  applicability  of  the 
perturbation  solution,  with  which  we  wish  to  compare  the  numerical  solution,  the  constant 
A  was  set  to  produce  a  relatively  low  peak  strain  amplitude  of  0.6  ^.strain.  Numerical 
solutions  were  obtained  by  both  finite  difference  and  pseudospectral  methods,  using  a  grid 
spacing  of  1/20  the  dominant  wavelength.  The  analytical  solution  for  this  source  was 
approximated  using  the  first  and  second  perturbation  terms  obtained  using  the  method 
given  by  McCall. 

We  specify  the  source  strength  in  terms  of  the  peak  strain  of  the  induced  plane 
wave.  We  present  results  in  dimensionless  form,  giving  numerical  results  for 
dimensionless  stress  S/Gq  and  velocity  u/c  (where  c  is  the  infinitesimal-strain  wavespeed). 
It  then  follows  from  the  scaling  relationship  (7)  that  the  solutions  for  these  dimensionless 
quantities  depend  on  propagation  distance  only  through  the  ratio  of  distance  to  wavelength 
of  the  dominant  frequency  (i.e.,  xlcT,  where  c  is  the  (infinitesimal-strain)  wavespeed),  and 
on  time  only  through  the  ratio  of  time  to  the  source  period  (i.e.,  t/T).  The  source  duration 
w  is  set  to  3  for  all  computations. 

Fig.  3  shows  the  velocity  spectra  at  distances  from  the  source  of  8  and  12  times  the 
wavelength  of  the  source  frequency.  The  left  half  of  the  figure  compares  the  numerical  and 
analytical  solutions  near  the  fundamental  frequency.  The  right  half  of  the  figure  shows  the 
comparison  of  solutions  at  the  second  multiple  (/=2/T)  and  third  multiple  (f=HT)  of  the 
source  frequency.  The  comparison  shows  that  the  two  numerical  methods  give  nearly 
identical  results  for  excitation  of  higher  harmonics,  and  that  both  agree  remarkably  well 
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with  the  perturbation  solution.  The  comparison  confirms  the  applicability  of  our  finite 
difference  and  pseudospectral  methods  to  nonlinear  problems,  and  incidentally  verifies  that 
McCall’s  perturbation  solution  is  highly  accurate  in  this  instance. 

We  perform  similar  computations  and  examine  the  evolution  with  distance  of  the 
harmonics,  and  their  dependence  on  source  amplitude.  For  these  computations,  the  source 
wave  shape  was  slightly  different  from  that  used  for  the  above  comparison;  we  used  a 
moment  density  to  excite  the  plane  wave,  resulting  in  a  velocity  waveform  which  is  the 
derivative  of  (28).  For  this  narrowband  source,  the  difference  is  negligible  in  both  time 
domain  waveform  and  in  the  spectral  envelope  (once  we  scale  the  source  amplitude  to 
produce  a  given  peak  strain  amplitude).  In  particular,  the  ratio  cited  above  between  time 
domain  peak  and  spectral  peak  (of  the  fundamental  harmonic)  still  applies.  The  stress 
waveforms  and  spectra  at  distances  4, 8, 12, 16  wavelengths  from  the  source,  respectively, 
are  shown  in  Fig.  4a  for  the  case  where  peak  strain  is  equal  to  1.2  p,strain.  As  expected, 
the  excitation  of  the  higher  harmonics  is  greater,  at  a  given  propagation  distance,  than  was 
the  case  for  the  lower  amplitude  case  in  Figure  3. 

Figures  4b  shows  the  variation  with  propagation  distance  of  both  velocity 
waveforms  and  spectral  amplitudes  for  the  case  of  the  1.2  pstrain  source.  Though 
significant  higher  order  harmonics  are  evident  in  the  spectra,  and  their  amplitudes  increase 
systematically  with  propagation  distance,  the  waveforms  themselves  are  nearly 
indistinguishable  at  the  different  distances.  The  inset  shows  the  amplitudes  of  the  second 
and  third  harmonics,  as  functions  of  propagation  distance,  in  the  form  of  ratios  to  the 
fundamental  (source  frequency)  harmonic.  Note  that  the  growth  rate  of  the  second  multiple 
approaches  an  asymptotic  slope  of  1,  i.e.,  proportional  to  distance.  This  can  be  compared 
with  the  approximate  theoretical  result  for  purely  monochromatic  excitation  (Thompson  and 
Tiersten,  1977;  McCall,  1994;  Meegan  et  al.,  1993),  according  to  which  the  displacement 
amplitude  Uiiof  the  harmonic  at  twice  the  source  frequency,  2/o,  depends  on  displacement 
amplitude  Ui  at  the  source  frequency /o,  the  wavenumber  tcQ  of  the  source  frequency,  the 
nonlinearity  parameter  j8,  and  distance  x, 

(29) 

This  result  can  also  be  written  in  dimensionless  form  as  an  expression  for  the  ratio  U2flUf, 
in  terms  of  the  ratio  L  of  distance  to  source  wavelength  and  the  strain  amplitude  at  the 
fundamental  frequency /o , 


(30) 

Uj  2^  ^ 
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The  asymptotic  slope  for  the  2/o  harmonic  in  Fig.  4b  indicates  growth  with  the  first  power 
of  distance,  in  agreement  with  (30).  The  growth  rate  differs  somewhat  from  (30),  as  we 
would  expect,  since  the  source,  while  quite  narrowband,  is  not  purely  monochromatic 
(note  also  that  Eqn  30  gives  the  ratio  of  displacement  amplitudes,  and  the  velocity  ratio 
U^fjUf  appropriate  for  comparison  with  Fig.  4b  is  a  factor  of  2  higher).  The  asymptotic 
slope  for  the  3/o  harmonic  is  2,  which  also  agrees  with  theory  (McCall,  1994),  according 
to  which  proportional  to  (jSEfL  )2^ 

Figure  4c  shows  the  variation  of  the  velocity  spectral  amplitudes  with  source 
amplitude,  at  a  fixed  distance  of  13  wavelengths  from  the  source.  The  calculations  are  again 
consistent  with  the  scaling  expected  on  the  basis  of  (30).  The  fundamental  component 
decreases  slightly  as  it  loses  energy  to  the  higher  harmonics.  The  ratio  U2flU f  increases 
nearly  linearly  with  source  strain  amplitude,  as  predicted  by  (30);  the  ratio  UyjUf 
increases  approximately  with  the  square  of  the  source  strain  amplitude,  also  in  agreement 
with  theory. 

Berea  Sandstone  Model  1— Strain-Independent  Elastic  Modulus 

The  nonlinear  elastic  simulations  in  the  last  section  were  useful  because  they  could 
be  compared  with  analytical  approximations,  and  that  comparison  demonstrates  that  the 
numerical  methods  can  be  used  successfully  to  model  pulse  propagation  and  harmonic 
distortion.  However,  the  nonlinear  elastic  model  is  not  realistic  for  rocks,  in  that  it  neglects 
hysteresis  and  attenuation.  In  this  section  we  examine  the  evolution  of  a  plane  wave 
propagating  through  the  Berea  sandstone  model  described  earlier,  represented  with  the 
endochronic  formulation.  This  section  considers  the  Berea  model  developed  from  quasi¬ 
static  data  alone,  which  has  a  constant  (strain-independent)  elastic  modulus.  Comparison 
of  the  results  to  data  from  wave  propagation  experiments  will  demonstrate  the  need  (in  the 
context  of  the  constitutive  formalism  used  here)  to  add  a  weakly  strain-dependent  term  to 
the  elastic  modulus. 

We  use  the  same  plane  wave  source  described  in  the  previous  section,  giving  a 
plane  wave  at  the  source  whose  velocity  waveform  is  the  derivative  of  (28).  Figure  5 
shows  the  stress  pulse  after  propagation  by  4,  8,  12  and  16  fundamental-frequency 
wavelengths,  for  the  case  of  a  source  with  peak  strain  of  1.2  pstrain.  While  little  change  in 
wave  shape  is  detectable  in  the  time  domain,  a  distinct  decay  of  peak  time-domain 
amplitude  is  evident,  in  contrast  to  the  nonlinear  elasticity  results  in  Fig.  4a.  Harmonic 
distortion  effects  are  also  evident,  as  seen  in  the  spectral  plots  in  Figure  5.  As  the  pulse 
propagates,  the  amplitude  of  the  fundamental  component  gradually  diminishes  with  the 
propagation  distance.  The  losses  in  the  fundamental  are  much  greater  than  the  very  small 
losses  (attributable  to  energy  transfer  to  higher  harmonics)  seen  in  the  fundamental  in  the 
nonlinear  elastic  examples  in  the  previous  section,  and  are  a  consequence  of  the  hysteresis 
in  the  Berea  model.  No  viscous  loss  mechanism  has  been  added  to  the  Berea  model  derived 
earlier  from  fitting  the  endochronic  model  to  laboratory  data.  Harmonics  representing  the 
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odd  multiples  of  the  fundamental  frequency  are  generated  and  progressively  increase  in 
amplitude  with  propagation  distance.  The  amplitudes  of  these  odd  harmonics  are 
somewhat  greater  than  they  were  for  the  nonlinear  elastic  case  at  the  same  strain  level  and 
propagation  distance  (i.e..  Fig  4a).  On  the  other  hand,  their  growth  rate  appears  to 
diminish  with  propagation  distance.  This  can  be  understood  as  a  consequence  of 
increasing  hysteretic  loss  at  higher  amplitude.  The  amplitudes  of  the  harmonics  are 
determined  by  competing  effects:  growth  due  to  energy  transfer  from  lower  harmonics  and 
attenuation  due  to  amplitude-dependent  losses  in  hysteresis. 

The  most  striking  characteristic  of  Fig.  5,  in  comparison  with  Fig.  4a,  is  the 
complete  absence  of  the  even  harmonics  in  the  Berea  simulations.  Spectra  from  wave 
propagation  experiments  in  Berea  sandstone  presented  by  Meegan  et  al  (1993),  do  show 
that  the  odd  harmonics  dominate  the  spectrum.  These  authors  propose  that  higher  order 
terms  (i.e.,  cubic  anharmonicity)  might  be  required  to  explain  the  relative  strength  with 
which  the  odd  harmonics  are  observed.  Our  simulations  show  that  incorporation  of  realist 
hysteresis  into  the  Berea  model  is  capable  of  explaining  dominance  of  the  odd  harmonics. 

On  the  other  hand,  the  spectra  of  Meegan  et  al.  still  show  significant  excitation  of 
the  even  harmonics,  even  though  the  odd  harmonics  are  dominant.  Other  observations 
conducted  in  the  intermediate  strain  regime,  from  both  laboratory  and  field  experiments, 
commonly  show  excitation  of  harmonics  of  both  even  and  odd  order  (e.g.,  Beresnev  and 
Nikolaev,  1988).  The  absence  of  the  even  harmonics  in  the  Berea  simulations  shows  that, 
despite  accounting  well  for  the  shape  of  laboratory  stress-strain  curves,  our  Berea  model 
lacks  the  slight  asymmetry  of  response  between  loading  and  unloading  directions  present 
in,  for  example,  the  nonlinear  elastic  model  with  cubic  anharmonicity.  This  model 
deficiency,  while  too  subtle  to  be  readily  evident  in  comparisons  of  quasistatic  stress-strain 
curves,  becomes  striking  evident  in  the  spectra  of  the  wave  propagation  experiments.  In 
this  respect,  the  generation  of  harmonics  provides  a  particularly  sensitive  test  of  nonlinear 
constitutive  models. 

Berea  Sandstone  Model  2-Strain-Dependent  Elastic  Modulus 

In  view  of  the  foregoing  results,  we  modify  the  Berea  model  through  the  addition 
of  a  strain  dependence  to  the  modulus  G,  as  in  (27),  again  using  /^5000,  and  keeping  the 
infinitesimal-strain  value  of  the  modulus.  Go,  equal  to  8.96x10^.  For  strains  up  to  about 
10‘5,  this  change  affects  the  modulus  by  at  most  5%.  The  modified  Berea  model  yields 
hysteresis  loops  whose  agreement  with  the  laboratory  data  is  indistinguishable  from  that  of 
the  model  with  strain-independent  modulus.  Yet  the  harmonics  generation  predicted  by  the 
model  is  affected  substantially. 

Figure  6a  shows  the  stress  pulse  at  distances  of  4,  8,  12, 16  wavelengths,  together 
with  the  corresponding  amplitude  spectra,  for  a  source  with  peak  strain  of  1.2  pstrain.  The 
spectra  exhibit  some  of  the  same  features  seen  with  the  preliminary  Berea  model,  including 
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the  progressive  attenuation  of  the  amplitude  of  the  fundamental,  strong  excitation  of  the  odd 
harmonics,  and  apparent  reduction  of  the  growth  rate  of  the  harmonics  with  propagation 
distance.  In  contrast  to  the  preliminary  model,  however,  the  modified  Berea  model  also 
excites  the  even  harmonics.  Figures  6b  shows  the  variation  of  particle  velocity  spectral 
amplitude  with  propagation  distance  (the  inset  shows  ratios  of  the  fundamental,  second, 
and  third  harmonics  to  the  initial  amplitude  of  the  fundamental).  The  second  harmonic 
increases  nearly  linearly  with  distance  and  the  third  harmonic  also  increases,  but  at  a  slower 
rate.  Comparison  of  Fig.  4b  with  the  corresponding  figure  for  nonlinear  elasticity.  Fig. 
4b,  shows  some  significant  differences.  The  Berea  model  excites  a  significantly  stronger 
3/o  harmonic,  for  example.  But,  as  a  result  of  hysteretic  losses,  the  amplitude  of  that 
harmonic  quickly  saturates  with  propagation  distance  in  the  Berea  model,  rather  than 
continuing  to  grow.  The  2/o  harmonic  also  shows  some  pronounced  effects  from 
hysteresis.  It  is  somewhat  weaker  at  all  distances  in  the  Berea  model,  compared  with  its 
amplitudes  in  the  nonlinear  elasticity  model,  even  though  the  elastic  modulus  as  a  function 
of  strain  is  the  same  in  both  cases.  The  second  harmonic  also  shows  a  diminishing  growth 
rate  with  increasing  distance  of  propagation. 

Both  second  and  third  harmonics  exhibit  approximately  the  same  rate  of  increase 
with  source  amplitude  in  the  modified  Berea  model,  as  Fig.  6c  shows.  Comparison  with 
Fig.  4c  again  shows  significant  differences  between  the  Berea  model  and  nonlinear 
elasticity.  The  amplitude  of  the  fundamental  is  nearly  proportional  to  source  amplitude  for 
the  nonlinear  elasticity  model,  but  diminishes  with  increasing  source  amplitude  in  the  Berea 
model.  The  latter  effect  is  a  result  of  nonlinearity  of  the  hysteretic  loss  mechanism,  which 
leads  to  amplitude-dependence  of  apparent  Q,  The  second  and  third  harmonics  have  quite 
different  rates  of  growth  with  source  amplitude  in  the  nonlinear  elasticity  case,  but  very 
similar  dependence  on  source  amplitude  in  the  Berea  model.  The  behavior  of  the  Berea 
model  reflects  the  competing  effects  of  amplitude-dependent  attenuation  and  energy  transfer 
among  harmonics. 

Experiments  have  shown  that  propagation  distortion  from  a  dual  frequency  source 
is  qualitatively  different  from  what  is  seen  with  a  single  frequency  source  (Beresnev  and 
Nikolaev,  1988;  Bonner  and  Wanamaker,  1990,  1991;  Johnson  et  al,  1987;  Johnson  and 
Shankland,  1989;  Johnson  et  al.,  1991).  While  a  single-frequency  source  incurs  only 
harmonic  distortion,  a  two-frequency  source  is  subject  not  only  to  harmonic  distortion  but 
also  to  intermodulation  distortion,  that  is  the  production  of  components  at  frequencies 
which  are  the  sums  and  differences  of  integer  multiples  of  the  two  source  frequencies. 

We  numerically  simulate  an  intermodulation  experiment  for  the  Berea  model  using  a 
two-frequency  source  consisting  of  a  sum  of  narrowband  sources  of  the  form  used  in  the 
previous  simulations.  The  two  source  frequencies  are  in  the  ratio  3:2;  their  displacement 
amplitudes  are  equal,  imparting  a  3:2  ratio  to  their  strain  (and  velocity  and  stress) 
amplitudes.  The  sources  are  scaled  to  give  a  combined  peak  strain  at  the  source  of  1.2 
pstrain.  Figure  7  shows  the  waveforms  and  the  corresponding  spectra,  at  distances  of  4, 
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8, 12,  and  16  wavelengths  from  the  source,  respectively.  Spectral  peaks  associated  with 
the  sum  and  difference  frequencies  (e.g.,  fi+f2, 2fi+f2,  fi-f2),  are  seen  in  this  picture,  as 
are  the  harmonics  of  the  two  source  frequencies  (e.g.,  2fi,  2f2).  At  the  relatively  low 
source  strain  amplitude  used  in  this  simulation,  the  harmonics  and  intermodulation  terms 
remain  relatively  small  compared  with  the  fundamental  source  frequencies.  Field 
experiments  have  shown  that  it  is  possible,  under  some  circumstances  and  for  some 
materials,  for  the  higher  harmonics  from  a  monochromatic  source,  or  the  interaction- 
generated  peaks  from  a  dual  frequency  source,  to  exceed  the  amplitudes  of  the  fundamental 
frequency  components  ( Beresnev  and  Nikolaev,  1988). 

3.7  DISCUSSION 

The  foregoing  model  for  nonlinear  rock  behavior  in  the  intermediate  strain  regime 
has  relatively  few  parameters:  Go,  P,  cc,  and  C.  Of  these,  the  first  is  simply  the  elastic 
modulus.  Of  the  3  nonlinearity  parameters,  P  is  just  the  conventional  cubic  anharmonicity 
parameter,  and  we  fix  a  to  the  value  1/2  to  match  the  generally  observed  power  law  scaling 
of  Qr^  with  strain  amplitude.  Only  C  was  adjusted  to  match  the  quasistatic  observations  of 
hysteresis  in  Berea  sandstone.  However,  C  can  also  be  related  to  the  slope  of  Q'^  versus 
strain  amplitude.  Thus,  the  model  may  be  applicable  in  cases  where  one  has  an  estimate  of 
this  slope,  either  from  laboratory  wave  propagation  experiments  or  from  seismic 
observations,  even  if  direct  measurements  of  stress-strain  curves  are  lacking.  That  is,  it 
may  have  utility  as  a  generic  model  for  nonlinear  seismic  wave  propagation  in  the 
intermediate  strain  range,  parameterized  by  the  anharmonicity  coefficient  and  the  Q'^  slope. 

Besides  having  parameters  which  can  be  related  to  fairly  standard  observables,  the 
model  has  the  additional  advantage  of  having  governing  equations  which  are  easily 
programmed,  require  minimal  storage,  and  are  efficiently  and  stably  solved.  These  latter 
attributes  are  essential  in  any  constitutive  model  which  is  to  be  used  in  large,  multi¬ 
dimensional  wave  propagation  simulations.  As  an  example,  a  typical  3D  simulation  with 
300^  grid  cells  will  require  of  the  order  of  10^®  to  10^^  evaluations  of  the  constitutive 
functional. 

The  endochronic  formulation  is  not  the  only  available  theoretical  framework  for 
modeling  the  phenomenology  of  plasticity  without  a  yield  surface  (e.g.,  Defalias  and 
Popov,  1976;  Bardet,  1996),  and  it  is  likely  that  other  formulations  can  also  be  adapted  to 
reproduce  key  observables  associated  with  intermediate  strain  levels  is  rock,  at  least  over 
some  limited  strain  range.  The  endochronic  formulation  may  itself  be  quite  limited  in  its 
strain  range  of  applicability.  Its  most  significant  limitation  is  probably  that  the  endochronic 
constitutive  model  inevitably  produces  so-called  "ratcheting"  when  cycled  repeatedly  over 
highly  asymmetrical  paths.  That  is,  high-amplitude  hysteresis  loops  cycled  repeatedly 
through,  e.g.,  one-sided  loads,  do  not  retrace  themselves  exactly.  This  phenomenon  of 
ratcheting  is  also  sometimes  called  cyclic  creep.  While  this  phenomenon  can  be  turned  to 
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advantage  for  modeling  the  response  of  granular  soils,  where  ratcheting  is  frequently 
observed  experimentally,  it  is  a  decided  limitation  for  modeling  consolidated  rock. 
Ratcheting  in  the  Berea  model,  for  example,  becomes  significant  at  strains  exceeding 
several  lO's  of  jistrain,  above  which  the  endochronic  formulation  may  not  be  appropriate 
for  some  applications. 

The  Fourier  pseudospectral  and  fourth-order  finite  difference  methods  work  well, 
without  the  incorporation  of  artificial  viscosity,  at  the  strain  levels  investigated  here. 
However,  this  is  largely  because  at  the  low  strain  levels  which  we  simulated,  there  was 
negligible  harmonic  generation  near  the  spatial  Nyquist  wavenumber.  At  higher  strains,  the 
stronger  nonlinearity  will  excite  those  much  higher  harmonics  significantly.  It  may  be  that 
incorporation  of  a  realistic  anelastic  mechanism  in  series  with  the  endochronic  model  will 
suffice  to  suppress  these  harmonics.  Otherwise,  it  will  be  necessary  to  use  smaller  grid 
elements,  artificial  viscosity,  or  both,  to  achieve  stable  simulations  at  higher  strain  levels.  If 
the  generation  of  higher  harmonics  is  not  sufficiently  counterbalanced  by  losses  in 
hysteresis,  it  may  become  desirable  to  apply  differencing  methods  tailored  to  strongly 
nonlinear  waves  (e.g.,  Cheng,  1996). 

Quasi-harmonic  numerical  simulations  with  the  Berea  model  reveal  wave 
propagation  phenomena  which  are  absent  in  the  nonlinear  elastic  simulations  and  theory. 
Decay  of  the  fundamental  frequency  component  was  very  evident  in  these  simulations,  in 
contrast  to  the  nonlinear  elastic  case.  In  the  nonlinear  elastic  simulations,  decay  of  the 
fundamental  was  negligible  at  the  strain  levels  and  propagation  distances  considered  here, 
resulting  only  from  energy  transfer  to  the  higher  harmonics.  In  the  Berea  model,  however, 
decay  of  the  fundamental  was  pronounced  and  strongly  amplitude  dependent.  The 
attenuation  of  the  fundamental  was  dominated  by  amplitude  dependent  hysteretic  losses 
rather  than  by  energy  transfer. 

Strong  excitation  of  the  third  harmonic  I/sf  was  also  evident  in  our  Berea 
simulations.  This  result  is  consistent  with  laboratory  observations  of  wave  propagation  in 
Berea  sandstone.  Johnson  and  Rasolofosaon  (1996)  suggested  that  cubic  anharmonicity 
alone  may  be  inadequate  to  explain  the  high  excitation  level  of  harmonic  U^t,  and  proposed 
that  hysteresis  effects  may  be  the  most  physically  reasonable  explanation  (higher  order 
anharmonicity  being  another  potential  contributor).  Our  Berea  simulations  provide  some 
support  for  this  proposal,  demonstrating  that  a  realistic  model  of  hysteresis  enhances  the 
relative  level  of  third  harmonic  excitation  in  a  manner  consistent  with  laboratory 
observations. 

The  Berea  simulations  also  exhibit  saturation  effects  in  the  growth  of  harmonics, 
i.e.,  reduction  of  harmonic  growth  rates  with  distance  of  propagation.  This  result  reflects 
the  amplitude  dependence  of  the  hysteretic  loss  mechanism-as  amplitudes  increase, 
increasing  attenuation  counteracts  the  energy  transfer  from  low  to  higher  harmonics. 
Similar  saturation  effects  for  harmonics  have  been  measured  experimentally  in  other  solids 
(e.g.,  in  copper  by  Breazeale  and  Ford,  1965). 
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It  is  easy  to  formally  generalized  the  endochronic  constitutive  model  to  2D  or  3D, 
and  the  numerical  procedures  will  not  be  significantly  more  complex.  Such  a 
generalization,  however,  entails  more  unknown  parameters,  including  parameters 
representing  possible  nonlinear  coupling  among  the  strain  components.  We  have 
successfully  fit  the  parameters  of  a  model  of  this  type  to  laboratory  observations  involving 
coupled  volumetric  and  shear  strains.  The  main  difficulty  is  that  there  is  no  guarantee,  a 
priori,  that  the  resulting  model  (or  any  other  purely  phenomenological  constitutive  model) 
will  accurately  predict  stress-strain  behavior  over  loading  paths  much  different  from  those 
used  to  obtain  the  fit.  This  caveat  will  apply  to  any  purely  phenomenological  constitutive 
model.  At  a  minimum,  the  method  continues  to  provide,  in  2D  and  3D,  a  convenient  way 
to  parameterize  laboratory  hysteresis  observations  in  a  form  suitable  for  efficient  wave 
propagation  computations.  It  may  be  possible  to  use  observations  of  wave  interactions 
(e.g.,  P-SV  wave  interactions,  as  in  Qian,  1995)  to  provide  further  constraints  on  2D  or  3D 
forms  of  the  model. 

3.8  CONCLUSIONS 

We  have  demonstrated  a  method  for  modeling  nonlinear  wave  propagation  in  rock 
that  appears  to  provide  a  stable,  efficient,  and  accurate  algorithm  for  numerical  simulations 
at  intermediate  strain  levels.  The  method  is  based  upon  a  singular-kernel  endochronic 
plasticity  model  of  hysteresis,  linked  to  either  a  finite  difference  or  pseudospectral  wave 
propagation  scheme.  The  constitutive  model  preserves  rate-independence  (and  simple 
scaling),  produces  hysteresis  loops  with  initial  loading  and  unloading  at  the  elastic  slope  at 
reversal  points,  and  provides  a  power  law  dependence  of  apparent  on  strain  amplitude. 
We  have  shown  that  the  model  parameters  can  be  successfully  fit  to  quasistatic  stress-strain 
curves  obtained  at  intermediate  strain  levels,  and  demonstrate  the  fit  for  Berea  sandstone. 
The  method  should  be  applicable  to  other  rocks  under  conditions  for  which  attenuation  is 
amplitude  dependent,  and  will  help  fill  an  existing  gap  in  seismic  theory  at  intermediate 
strains. 

Application  of  the  Berea  sandstone  model  to  quasi-harmonic  wave  propagation 
reveals  a  number  of  contrasts  with  comparable  simulations  using  nonlinear  elasticity  (with 
cubic  anharmonicity).  (i)  The  amplitude  of  the  fundamental  frequency  component  decays 
much  more  rapidly  in  the  Berea  model,  indicating  significant  amplitude-dependent, 
hysteretic  losses,  (ii)  The  harmonic  at  3  times  the  fundamental  frequency  is  significantly 
larger  for  the  Berea  model,  especially  at  short  propagation  distances,  and  in  better 
agreement  with  laboratory  measurements,  (iii)  Growth  of  higher  harmonics  with 
propagation  distance  in  the  Berea  model  shows  a  saturation  effect,  a  behavior  which  is 
consistent  with  amplitude  dependent  attenuation,  but  which  is  absent  from  the  nonlinear 
elasticity  model.  Application  to  dual-frequency  wave  propagation  predicts  the  excitation  of 
intermodulation  distortions,  with  the  production  of  components  corresponding  to  sums  and 
differences  of  integer  multiples  of  the  source  frequency. 
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The  methods  of  this  paper  have  potential  application  to  both  applied  and 
fundamental  research  problems  in  seismology,  including  modeling  the  influence  of  near¬ 
source  intermediate  strains  on  the  source  signature  of  earthquakes  and  underground 
explosions.  The  methods  may  also  have  applications  to  the  development  of  nondestructive 
seismic  methods  for  in  situ  rock  mass  characterization. 
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FIGURE  CAPTIONS 


Figure  1:  Fitting  the  constitutive  model  to  Berea  sandstone  data,  (a)  Calculation  of  the 
elastic  modulus  G  and  hysteresis  parameter  C  from  the  reversal  points.  Stars  denote 
observations,  and  the  dashed  curve  is  the  model  prediction  using  the  optimal  values  of  G 
and  C.  The  parameter  a  was  fixed  at  a  value  of  l/2.(b)  Prony  approximation  to  the  kernel 
function,  in  the  range  z=2xl0-*  to  4x10-^.  The  exact  kernel  is  shown  as  a  solid  curve,  the 
approximation  as  a  dashed  curve,  (c)  The  resulting  synthetic  hysteresis  loops  calculated  by 
direct  integration  using  the  exact  kernel  function,  at  maximum  strain  levels  of  10, 15, 20, 
25,  30,  35, 40, 45  and  50  jxstrains,  respectively,  (d)  Same  as  (c),  but  calculated  using  the 
Prony  approximation  to  the  kernel. 

Figure  2:  Comparison  of  hysteresis  loops  calculated  from  the  Prony  approximation 
kernel  described  in  Fig.  1,  dashed  lines),  with  observed  loops  for  Berea  sandstone,  solid 
lines  (G.  Boitnott  and  R.  Haupt,  New  England  Research,  Inc.,  personal  communication). 
Peak  strain  levels  are  7.4, 9, 19.6  and  37  ^strains,  respectively. 

Figure  3:  Propagation  of  a  quasi-harmonic  plane  wave  in  a  nonlinearly  elastic  medium: 
comparison  of  approximate  analytic  solution,  obtained  by  the  method  of  McCall  (1994), 
with  finite  difference  and  pseudospectral  numerical  solutions.  The  analytic  solution  is 
shown  by  the  solid  curve.  The  two  numerical  solutions  are  indistinguishable  from  each 
other,  and  only  the  pseudospectral  case  is  plotted.  The  source  was  scaled  to  give  a 
maximum  time  domain  strain  of  0.6x10"^.  Diamonds  and  stars  denote  the  numerical 
solution  at  distances  of  8  and  12  times  the  source  wavelengths,  respectively.  The  left 
portion  of  the  figure  shows  the  spectrum  near  the  source  frequency,  and  the  right  portion 
shows  the  next  two  higher  order  harmonics.  Frequency  has  been  normalized  by  the  source 
frequency  (1/7).  The  velocity  spectral  amplitude  V  has  been  normalized  by  the  wavelength 
of  the  source  excitation  (l/cT). 

Figure  4:  (a)  Numerical  solution  for  the  evolution  of  a  modulated  sinusoidal  pulse,  with 
increasing  propagation  distance,  in  a  nonlinearly  elastic  solid  with  cubic  anharmonicity. 
The  maximum  strain  of  the  source  pulse  is  1.2x10"^.  Stress  waveforms  (shown  in  the 
inset)  and  the  corresponding  spectra  are  shown  for  propagation  distances  of  4,  8,  12,  and 
16  wavelengths,  (b)  The  same  results  shown  in  the  form  of  particle  velocity  spectra.  The 
inset  shows  the  spectral  ratios  of  the  2/o  and  3/o  harmonics,  respectively,  to  the 
fundamental.  The  harmonics  increase  with  distance  as  predicted  by  theory,  (c)  Variations 
of  the  amplitude  spectra  with  strain  amplitude  of  the  source,  at  fixed  position  13 
wavelengths  from  the  source.  The  amplitude  scaling  of  the  numerical  solution  agrees  with 
theory.  The  source  strain  amplitudes  are  0.24, 0.6, 1.2,  and  2.3  pstrain,  respectively. 
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Figure  5:  Numerical  solution  for  the  evolution  of  a  modulated  sinusoidal  pulse,  with 
increasing  propagation  distance,  in  the  preliminary  Berea  model  (i.e.,  without  strain 
dependence  of  the  elastic  modulus).  The  maximum  strain  of  the  source  pulse  was  1.2xlOA 
Waveforms  (shown  in  the  inset)  and  the  corresponding  spectra  are  shown  for  propagation 
distances  of  4,  8,  12,  and  16  wavelengths.  Note  the  complete  absence  of  even-order 
harmonics  in  this  case,  as  well  as  the  attenuation  of  the  fundamental  component  with 
distance,  distinguishing  this  model  from  the  nonlinear  elastic  model. 

Figure  6:  (a)  Numerical  solution  for  the  evolution  of  a  modulated  sinusoidal  pulse,  with 
increasing  propagation  distance,  in  the  final  Berea  sandstone  model  (with  strain-dependent 
elastic  modulus).  The  maximum  strain  of  the  source  pulse  is  1.2x10'^.  Stress  waveforms 
(shown  in  the  inset)  and  the  corresponding  spectra  are  shown  for  propagation  distances  of 
4,  8, 12,  and  16  wavelengths,  (b)  The  same  results  shown  in  the  form  of  particle  velocity 
spectra.  The  inset  shows  the  spectral  ratios  of  the  2/o  and  3/o  harmonics,  respectively,  to 
the  fundamental.  The  increase  with  distance  of  the  higher  order  harmonics  shows 
saturation  effect  associated  with  nonlinear  attenuation,  (c)  Variations  of  the  amplitude 
spectra  with  strain  amplitude  of  the  source,  at  fixed  position  13  wavelengths  from  the 
source.  The  source  strain  amplitudes  are  0.3, 0.6, 1.2,  and  2.4  |i.strain,  respectively. 

Figure  7;  Numerical  solution  for  the  evolution  of  a  sum  of  modulated  sinusoids,  with 
increasing  propagation  distance,  in  the  final  Berea  sandstone  model  (with  strain-dependent 
elastic  modulus).  The  maximum  strain  of  the  source  pulse  is  1.2xl0'6.  Waveforms 
(shown  in  the  inset)  and  the  corresponding  spectra  are  shown  for  propagation  distances  of 
4, 8, 12,  and  16  wavelengths.  Note  the  generation  of  spectral  peaks  a  frequencies  equal  to 
sums  and  differences  of  integer  multiples  of  the  source  frequencies. 
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Table  1.  Model  Parameters 


Constitutive  Model 


Go 

8.96  X  109  Pa 

C 

3.87  X  107  Pa 

P 

5x  103 

a 

1/2 

Prony  Expansion 


Ai 

3.61  X  lOlO  Pa 

A2 

1.49  X  lOll  Pa 

A3 

5.67  X  1010  Pa 

A4 

4.56x1011  Pa 

ai 

1.0  x  105 

ai 

5.07  X  106 

^3 

2.75  x  107 

04 

9.27  X  107 
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FIGURE  1  (b) 
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FIGURE  2 

FIGURE  4  (a) 
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Normalized  Velocity  Spectrum  (V/cT) 


Normalized  Frequency  (fT) 


FIGURE  4  (b) 
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FIGURE  5 


Normalized  Stress  Spectrum  (S/GT) 


Normalized  Frequency  (fT) 


FIGURE  6  (a) 
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Normalized  Velocity  Spectrum  (V/cT) 


Normalized  Frequency  (fT) 


FIGURE  6  (b) 
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Normalized  Velocity  Spectrum(V/cT) 
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Normalized  Frequency  (fT) 


FIGURE  6  (c) 
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Normalized  Stress  Spectrum  (S/GT) 
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Normalized  Frequency  (fT) 


FIGURE  7 
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